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ABSTRACT 


A computer coded Lagrangian marker particle in Eulerian 
finite difference cell solution to the three dimensional 
incompressible mass transport equation. Water Advective 
Particle in Cell Technique, WAPIC, was developed, verified 
against analytic solutions, and subsequently applied in the 
prediction of long term transport of a suspended sediment 
cloud resulting from an instantaneous dredge spoil release. 
Numerical results from WAPIC were verified against analytic 
solutions to the three dimensional incompressible mass 
transport equation for turbulent diffusion and advection of 
Gaussian dye releases in unbounded uniform and uniformly 
sheared uni-directional flow, and for steady-uniform plug 
channel flow. WAPIC was utilized to simulate an analytic 
solution for non-equilibrium sediment dropout from an 
initially vertically uniform particle distribution in one 
dimensional turbulent channel flow. Predicted sediment 
fallout rates and vertical concentration values, as a 
function of downstream distance were within an absolute 
accuracy of four percent of the analytic values. 

The coastal zone mass transport version of the WAPIC 
algorithm employs a finite difference grid network which 
horizontally translates with the mean particle motion and 
expands with the dispersive growth of the marker particle 
cloud. The cartesian vertical co-ordinate of the three 
dimensional mass transport equation was transformed into 
dimensionless space, and then rearranged into flux 
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conseryative, pseudo total velocity form, in order to allow 
adaptation to flow situations with a time and spatiairy 
varying bottom topography and free surface. k velocity 
processor was developed for the interpolation of velocities 
from a vertically averaged two| dimensional velocity field 
of a tidally driven finite difference hydrodynamic 
numerical model. The transformed version of WAPIC was then 
coupled to the hydrodynamic velocity field by the velocity 
processor, and utilized to predict the long term diffusion 
and advection of dilute neutrally and negatively buoyant 
suspended sediment clouds resulting from an instantaneous 
release of barged dredge spoil at Brown’s Ledge in Rhode 


Island Sound 
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I. INTRODUCTION & STATEMENT OF PROBLEM 

Presently studies are in progress at the University of 
Rhode Island under the sponsorship of the U * S. Army Corps 
of Engineers concerning the feasibility and ecological 
impacts of the dumping of barged wastes in coastal areas 
such as Brown's Ledge. (1,2,3) From 1967 to 1970 
approximately 9.2 million cubic yards of dredge spoil was 
deposited within a one square nautical mile site, four 
miles south of Newport, in Rhode Island Sound. The work of 
Saila, Pratt» and Polgar states that, within the. limits of 
bathymetry measurements, practically little or no spoil was 
suspended during or after the disposal operation (2) . 
Presently there is concern that sediment suspended in the 
water column by vertical turbulence during and after the 
disposal operation may drift into Rhode Island Sound 
recreation areas. As a coastal zone management tool, the 
following work attempts to numerically model the long term 
transport (Ihr. - 2 day) of such wastes in the water 
column. WAPIC, a Lagranyian marker particle in Eulerian 
cell technique has been developed for simulating dilute 
sediment and particulate transport in a tidally dominated 
coastal zone environment. 

A literature review by Johnson (4) of the U. S. Army 
Corps of Engineers Waterways Experiment station, states - 
"Very little mathematical modeling of the physical fate of 
dredged material disposed of in an aquatic environment has 
been undertaken." Pertinent to this study are the 
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mathematical models of Koh and Chang (5) and Krishnappan 
(6), which consider the short term fate of the waste 
release cloud as it drops from the surface until it 
collects on the bottom. Disposal of dredged sediment or 
bottom material in the coastal zone is accomplished with 
barges. While dredged waste is normally released by an 
instantaneous hopper release, sewage sludges are pumped in 
a liquid slurry through a pipe while the barge is under 
way, and acids and other chemicals are released 
continuously in the wake of a moving barge. Koh and Chang 
produced a numerical model designed to simulate all three 
of the previously mentioned modes of discharge. Racge 
disposal of dredged sediment in Rhode Island Sound was done 
by an instantaneous dumping in relatively shallow coastal 
areas of depths between 90-110 feet (1) . Of interest is 
Koh and Chang's treatment of the instantaneous dumping 
operation . 

Their model assumes that the descent of the waste 
cloud takes place in three stages, each stage being 
governed by a different mechanism. They are- 

1. Convective cloud descent and entrainment 

2. Cloud collapse 

1. long term cloud transport settling ( Figure 1-1) 
During the first stage the negatively buoyant cloud 
descends and entrains water until neutral buoyancy is 
achieved; The assumption is made that the sediment laden 
cloud can be treated as a liquid of eguivalent density 
falling through a less dense ambient medium. The sediment 
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cloud which is now neutrally buoyant spreads horizontally 
at a constant vertical position due to liquid entrainment, 
in what is called the collapse stage. The assumption is 
made that the governing mechanism is the difference in 
hydrostatic pressure between the cloud and the ambient 
fluid. After the horizontal spreading stops the cloud of 
sediment particles undergoes vertical settling, and long 
term horizontal advection by the ambient current and 


turbulent diffusion 

due 

to this 

field . 

Koh and 

Chang * s 

model assumes that 

the 

ambient 

current 

field 

is 


horizontally unidirectional and steady, and may vary in 
magnitude only in the vertical direction. 

Krishnappan (6) developed a model for the determination 
of the size and height of a mound of sediment from an 
instantaneous release of dredge spoil in deep water with a 
turbulent one dimensional steady flow field. Laboratory 
experiments in still water showed that, the motion of 
uni formly sized particles can be treated in two stages, the 
initial entrainment phase, a rid the final settling phase. 
During the entrainment phase, the size of the particle 
cloud grows due to entrainment of ambient liquid While 
descent velocity decreases. During the settling phase, the 
fall velocity of the cloud is the same as the average 
settling velocity of the individual particles inside the 
cloud. In this stage, the cloud grows horizontally due to 
ambient turbulence, until it lands on the bottom. 

Krishnappan disagrees with the assumptions used by Koh 
and Chang in the first stage, and also .stated that the 
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collapse stage was unnecessary. Laboratory experiments 
showed that, during the entrainment or convective descent 
stage, the waste cloud must not be treated as a higher 
density liquid moving in an ambient fluid, but as a 
particle-liquid cloud or slurry moving through the water 
column. Results showed that the particle-liquid cloud did 
not increase in size throughout its first stage, while a 
liquid cloud of equivalent density continuously increases 
in horizontal size. The former grows until it reaches a 
constant size after which the cloud settles with decreasing 
velocity until it moves with the settling velocity of its 
particles. ( Figure 1-2) His results further demonstrate 
that, when the depth of the water is less than a thousand 
meters, the settling phase will not occur. Thus in the 
Great Lakes or in shallow estuarine or coastal zone waters, 
the possiblities that the settling phase may occur are 
rare, and the entrainment phase will dominate until the 
particle cloud lands on the bottom. Fie states - 
” Consequently models sucii as Koh and Chang's which devote 
considerable attention to the long term transport phase are 
unnecessarily elaborate." 

The physical insight derived from Krishanppan's work 
is that Isarged sedimeiitary waste that is dumped 
instantaneously in relatively shallow water, such as Blvode 
Island Sound waters, fails to the bottom in the 
"entrainment phase" with very little if any of the spoil 
remaining in the water column within a few hours of the 
disposal operation. Observations by Saila et. al. (1) 
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support this conclusion. 

Indeed a very conservative enqineering estimate, in 
Rhode Island Sound waters, can be made assuming that a 
maximum of 1-5 % of the total spoil volume remains 

entrained in dilute form in the water column immediately 
following the disposal operation by the action of ambient 
turbulence. This cloud of sedimentary particles possessing 
near neutral buoyancy with still water settling rates of 
0.001 to 0.03 ft/sec (for particle diameters ranging from 
0.00005 to 0.0002B ft. ) may remain in a water column 60 to 
150 ft. deep for a time span of one to twenty four hours. 
Should the ambient turbulence be great enougn to overcome 
the particle settling rates, the particles may be treated 
as if they were neutrally buoyant parcels of fluid. This 
cloud of visible turbidity will convect and diffuse with 
the ambient tidal and wind drift currents and turbulence to 
background levels. 

The question might be asked - Why develop a long term 
transport model for this 1-5 X of the dredged waste which 
may remain in the water column ? Isn't this amount of 
material too miniscule to demand simulation ? The answer 
is no. First of all, this single barge disposal is 
repeated at least a hundred times to deposit three million 
cubic yards of sediment. Thus 1-5 barge loads of sediment 
waste will be entrained by the ambient current field, and 
adyected and diffused across the coastal shelf, during this 
period of disposal operations. These dredged bottom 
sediments contain entrained pollutants in the form of 


8 . 


anaerobic sediments, heavy metals, PCB'S, etc. . Thus a 
biological oceanographer might be interested in the 
dispersal pattern of this cloud across the shelf as the 
particles eventually settle to the bottom. A single long 
term sediment transport simulation may provide estimates of 
the fallout patterns of this sediment and the entrained 
pollutants, ie., to demonstrate over what area the benthic 
environment is contaminated by the disposal operation. 

A prediction of fallout or dispersal patterns of 
dilute suspended sedimentary particles in a turbulent 
tidally driven advective field requires the use of a three 
dimensional mass transport model of WAprC's capability. 
What is required is the ability to model the dispersal and 
settling patterns of non-buoyant dilute particulate matter 
in a coastal or estuarine environment. 

A shortcoming of both Koh and Chang's and 
Krishnappan's work is the lack of proper treatment of the 
long term mass transport of the suspended plume by the 
spatial and temporally varying mean turbulent velocity 
field. In tidally driven waters such as Rhode Island 
vSound, the horizontal ambient velocity field is a 
complicated function of space and time. A cloud of 
suspended sediment particles will drift and diffuse with 
the tidal currents into areas where the flow field varies 
in magnitude and direction, from that found at the initial 
dump site. A long term transport treatment of suspended 
sediment in turbulent tidally dominated environment 
requires the coupling of a hydrodynamic velocity field with 
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a mass transport nodal. Spaulding and Hurlburt (8) and 
Hunter and Spaulding (7) have successfully coupled a two 
dimensional vertically averaged hydrodynamic model to two 
and three dimensional mass transport models, respectively^ 
for predicting transport of neutrally buoyant dissolved 
constituents in Narragansett Bay and Rhode Island Sound 
waters respectively. Since these techniques are strictly 
Eulerian mass transport models, these methods cannot 
describe the settling or deposition patterns of discrete 
particles with varying settling rates. The following work 
is directed towards developing a numerical model f!oc long 
term simulation (Ihr. to 2 day) of sedimentary transport 
and deposition in estuarine and coastal waters. No effort 
has been made to improve upon Koh and Chang's and 
Krishn appan' s 'short term' treatment of the descent of the 
dense sediment cloud to the bottom. 



WAPIC, a 

three dimensional 

Lagrangian marker 

particle 

in 

Eulerian 

cell technique has 

been 

developed to 

simulate 
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of dilute sedimentary 
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estuarine environment. 
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II. HISTORY of WAPIC 

During the last ten years a considerable amount of 
work has been done in the area of mathematical modeling of 
traiisport phenomena in coastal and estuarine waters (9, 10) . 
Usually these numerical models make the assumption that the 
particulate or pollutant mass is dilute within the flow 
field. Numerical mass transport models for estuarine or 
coastal zone waters are concerned with water quality 
parameters such as dissolved oxygen ( DO) , biochemical 
oxygen demand ( BOD), coliform decay, and larval growth 
rates and transport. (7,11,12) 

Models have been developed that solve the 
incompressible mass transport equation in either one, two, 
or three dimensions. The solution viewpoint may be either 
purely Eulerian, Lagrangian, or a combination Eulerian- 
Lagrangian. In the ideal Lagrangian approach, the fluid or 
flow field is assumed to be divided into numerous finite 
size zones each of which represents a parcel of fluid. 
This mesh of cells flows along as the fluid flows in a 
laanner characterized by the approximation to the partial 
differential equation of mass transport ( Harlow (13)). 
Pure Lagranqian models for estuarine or tidal pollutant 
transport have been developed by Fischer ( 1^) , and Wallis 
(15). These models seem to enjoy success when a 
complicated flow field may be broken up into one 
dimensional segments or volumes. 

The pure Eulerian scheme, divides the fluid field into 
a fixed mesh of cells or elements. This mesh is fixed in 


space, and th© fluid and particulates rapve throuc^ii 

/ o 

mesh., The Ruletian schemes may be subdiyided iit'tvO 
types of solution procedures, the finite ana 

finite element methods. The finite differeijce schemes 
divide the regime into rectangular grids with the dependent 
variables defined at either the cell corners, centers or 
faces. The partial differentials are linearly approximated 
across tiiese grids (16). Finite difference schemes may be 
further subdivided into the explicit and implicit methods 
of solution advancement in time. The explicit solution 
technique solves algebraically for the new value of the 
dependent variable from the previous value (16) t The 
implicit techniques solves for the new value a!t time t +At 
using advanced time values, requiring the use of iterative 
or Gaussian elimination techniques (16). Recently, in the 
field of estuarine pollutant transport modeling, two 
dimensional models have been popular, where the third 
dimensional independent variable has been comoved by 
integration techniques (11,17,18), The alternating 
direction implicit technique, ADI, has also enjoyed success 

r;; — 

in two dimensional hydrodynaaic and mass transport models 
by Leendertse (1^), Hess and White (18) , and Spaulding 
(20). The ADI method has been utilized by Hunter and 
Spaulding (7) in a three dimensional model for coliforra 
decay and transport in the Providence River . 

The second type' of Eulerian scheme, the finite element 
technique, seems particulacly well suited for modeling mass 
transport in situations with complicated boundaries that 
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fin-lte difference techniques might find difficult to 
resolve. This method divides the zone, of flow into 
triangular or quadrilateral elementsments with the 
depend^ht variables defined at the nodes of the triangles. 
A popular finite Sleaent technique solution technique is 
the Galer.kin method of weighted residuals. A review ;aiid 
explanation of this solution technique may be found in the 
works of Hang (21) , and leimkuhler (22) . Leimkuhler has 
recently compl^tsd ,a two dimensional vertically averaged 
finite element model for simulating sediment transport in 
Massachusetts Bay. ‘Kf one a Ariathurai (23) developed a 
One dimensional finite relement model for determining 
sediment transport irithe;Savannah River estuary. 

Interestingly, both finite element and finite 
difference Eulerian techniques share the following 
shortcomings in their computational procedures: 

1. Fictitious diffusion 

Fictitious diffusion is especially noticeable since each 
cell is forced to be homogeneous. When a material enters a 
cell all jits properties are uniformly diffused thrbuqhout 
the Eulerian cell. The B'agnitude of this, fictitious 
diffusion depends on the length and time scale of the 
numerical simulation. - 

2. Production of severly depressed concentration 
qrad^ients upstream of point source waste loads. i j 

3. Difficulty in handling scale dependent diffusion 
processes which are common in ocean and coastal zone 


regions. 
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4. Difficulty in representing an effluent of various 
particle sizes and density or buoyancy characteristics. 

The steps leading to a practical solution to these 
problems which are inherent in pure Eulerian schemes for 
turbulent mass transport can be traced in the development 
of Lagrangian- Eulerian numerical methods. 

In 1957, Evans and Harlow (24) of the Los Alamos 
Scientific Laboratory introduced the Particle in Cell, PIC, 
method. This numerical procedure was developed for the 
simulation of problems involving the dynamics of 
compressible fluids in two dimensions ( 13 ) . This hybrid 
technique requires two computational meshes or systems one 
Eulerian, the other Lagrangian. The domain through which 
the fluid flows is a fixed Eulerian mesh of computational 
cells. The fluid or mass is represented by particles with 
physical mass, energy, and velocity (etc.) which move 
through the Eulerian mesh. As stated by Roache (16) 
concerning the PIC method - The most unique aspect is 
that continuum flow is not modeled, rather a finite number 
of particles are used, their locations and velocities being 
traced by Lagrangian kinematics. They are not merely 
marker particles, but they participate in the calculation 
even when free surfaces and interfaces are present.*’ 

In 1963 the PIC method still had some shortcomings 
which Had not been overcome (13) . The PIC numerical method 
shared the Eulerian problem of resolving the detail of 
minute processes within the flow field, and a fixed grid 
system which could not overcome boundary difficulties. The 
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nost severe limitation of the method that increased the 
difficulty of expanding the method to simulate three 
diaensional flow processes was the large computational 
storage requirements and run times required for the 
Lagrangian and Eulerian meshes. The development of larger 
and faster computational systems should allow the 
adaptation of the method to more complicated flow 
simulations. 

The PIC method was strictly a solution technique for 
fluid dynamicists until sklarew (25) , in 1970, presented a 
paper detailing what shall be called the "pseudo total flux 
or velocity" method of solution for the turbulent mass 
transport equation* Sklarew's method allowed the 
adaptation of PIC for methods of solution to the turbulent 
three dimensional incompressible mass transport equation. 
The equation may be written in vector form as- 

^ + y.Uc - ^ . (K^jC) = 0 (2.1) 

and may be rewritten as 

|£+V-[(7-K..^)o]. 0 (2.2) 

or 

V* [UjCl = 0 {2.3) 
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where the pseudo total velocity vector is defined as 



= Diffusion velocity vector (2.4) 

U = Advective velocitjr vector 

where UL,. is the sunmation of the advective and diffusive 
velocity vectors, the the diffusion velocity being defined 
as , 

-K. . 

u = — ii T7 c (2.5) 

D c ^ 


Given an eKternally supplied velocity vector U, each 
Lagrangian Barker particle is moved explicitly in time with 
the total pseudo velocity, which is the summation of the 
advective and negative diffusion velocities. The diffusion 
velocity is calculated by the product of the marker 

f ■ 

particle concentration gradient across an Eulerian grid 
cell, divided by the concentraton at that positon, 
multiplied by the turbulent diffusion coefficient { Eg. 

2.5 ) . 

In 1971 Hotchkiss (26) applied an explicit finite 
difference solution technigue to the three dimensional 
incompressible Navier - Stokes equations to produce a time 
dependent advective atmospheric velocity field, Osing 
Lagrangian marker particles and Sklarew's pseudo total 
velocity technique, an atmospheric heat and mass transport 
model was developed, Sklarew and Hotchkiss recognized, as 
did Hirt and Cook (27) , that the treatment of minute detail 
of the sub flow field missed in the pure Eulerian schemes 



could now be overcome. Gravity, viscous drag as a function 
|j of particulate diameter, direct particle deposition, 

! reflection and scavenging from the boundaries could now be 

f applied with the use of Lagranqian marker particles for 

; simulation of airborne transport of pollutant particulates. 

i ' 

Recognizing the potential, Lange (28,29), Knox, Hardy, 

f '• 

t and Sherman (30) of the Lawrence Li verm ere Laboratory, 

i ■ 

coined a new phrase for this hybrid Eulerian - Lagranqian 
mass transport model, ADPlC. The Atmospheric Diffusion 
Particle in Cell technique was developed for prediction of 
atmospheric pollutant transport under complicated advective 
wind conditions, surface boundaries, and turbulence scales. 
This version of the PIC mass transport algorithm introduced 
an expanding grid system, borrowed from ALE, Alternating 
; Laqrangian Eulerian technique, developed by Hitt et. al. 

(31) for simulating dynamics of high speed flow with 
complicated boundary conditons. This grid system could 
translate and expand with the marker particle cloud as it 
' advected through space thus allowing greater resolution of 

particle cloud movement in space, and subsequently 
minimization of simulation computer time usage, since the 
simulation time step can be maximized. Scale dependent 
diffusion became a readily applicable feature by extracting 
a diffusion coefficient at each time step hy use of the 4/3 
power law equation which itself uses the horizontal 
standard deviations of the marker particle cloud. This 
feature will be explained in a subsequent chapter. 

Knox et. al. (30) proposed the utilization of the 
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A DP IC technique for the physical fate prediction of dredge 
spoil and application to modeling of three dimensional 
transport of sediment in San Fransisco Bay. After an 
extensive review of the literature of numerical modeling of 
dredge spoil disposal, Johnson (4) recommended the use of 
the ADPIC technique for simulation of coastal zone and 
estuarine sediment and pollutant transport. 

Spaulding (32) recognized the need for a three 
dimensional numerical model that was capable of handling 
the long term transport of discrete sedimentary particles 
in a coastal zone tidally dominated environment. WAPiC, 
Water Advective particle in Cell Technique, a three 
dimensional explicit finite difference, pseudo total 
velocity solution to the turbulent mass transport equation 
was developed by Spaulding (32) in 1975 at the NASA Langley 
Research Center. The following work is a continuation of 
the adaptation of Harlow and Evan Vs PIC method and 
Sklarew's pseudo total velocity technique for solving mass 
transport problems in the coastal zone environment. 
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III. DEVELOPMENT OE EQUATIONS 


A, Turbulent Advection Diffusion Equation 

The incompressible ensemble time averaged three 

! 

dimensional mass transport equation may be written as 




3y az jx ' xax ay y^ 


* ^ 5 ^) 

3z z az 


(3.1) 


u,v,w 


= x,y,z time and space averaged velocity components 


K ,K ,K = x,y,z turbulent diffusion coefficients 
X y z 


= Time and spatially averaged concentration 


This equation assumes that the pollutant is dilute and 
neutrally buoyant and that the flux terras due to molecular 
diffusivity can be neglected since, in riverine, oceanic and 
atmospheric modeling of mass transport, these terms are 
typically two to three orders of magnitude smaller than the 
turbulent coefficients (17,33,35) . For purposes of this 
work, the above equation also assumes that the fluid flow 
field may be specified by some time and spatially averaged 
velocity components and that the coefficients are time and 
spatially averaged "dispersion” coefficients which may be 
empirically or theoretically determined. 


B, Addition of Particle Settling Velocity 
In order to simulate the motion of a negatively 
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buoyant dilute constituent such as suspended sediment 
particles, Eq, 3.1 may be modified by the addition of the 
tern 



which is the advective mass flux in the Z-direction due to 
the still water asymptotic particle settling velocity, ffs. 
The assumptions are made that the free fall particle 
settling velocity is constant in time and space for a given 
particle diameter and that the suspended constituent is 
dilute within the flow field. Equation 3.1 can now be 
rewritten 


a|£ + v|°- + (k i^) 

X dx 


(K 

y 


4 (K^) -w 

ZQZ s ^Z 


(3.2) 


This equation, the turbulent incompressible three 
dimensional mass transport equation for negatively buoyant 
dilute particles is the basic equation around which the 
WftPTC solution method is constructed. 


C= Flux Conservative Form 

WAPIC employs the pseudo total velocity method which 
may be conceptualized as a flux conservative formulation o£ 
the three dimensional mass transport equation. Equation 

^ ^ BBTODUGEBIUTY- (>F THE 
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3,2 may be rearranged - 


4c . i '''x4o,o 4 4c o 


K 


^ 'z dc c 

+ ^ (W + W = 0 

^z sc ^z 


(3.3) 


and rewritten as 


+ j 3L m c) + -^ (V c) + ^ (w c) 
Jt 9x ' T ^ ' T ' Jz ^ T ^ 


= 0 


(3.4) 


where 


u„ = u + 

T D 


(3.5a) 


v„ = V + 

T D 


W = W 4- 
T 3 D 


(3.5b) 

(3.5c) 


and 


K 


= 


_X ^1 
c 


(3 .6a) 


K 


= 


x\- 

c~ fy 


(3.6b) 


K 


= 


z )c 

C 9z 


(3.6c) 


where U^, V^, are defined as the pseudo total transport 

velocities. Equation 3.4 states that a particle is 
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advQctPd with the mean motion o£ the fluid field and 
diffused with a velocity proportional to the concentration 
qratlient surrounding that particle. 


D. Transformation of the Vertical Coordinate 

The numerical solution to Eg. 3.4 requires a uniform 
mesn of equivalent shape and volume finite difference cells 
throughout the three dimensional cartesian grid system. 
When this system of equal volume cells must be applied to 
two and three dimensional hydrodynamic velocity fields with 
time and spatially varying free surfaces and irregular 
bottom topography complications arise. In estuaries and 
coastal areas, the bottom topography may be irregular and 
will not fit a single rectangular mesh of equal volume 
cells, especially when that mesh of equal volume cells is 
translating and expanding in time beneath a free surface 
which is also rising and falling with space and time. 

Following the method of Hunter and Spaulding (7) , and 
Gordon (36), a solution is to non-diraensionalize the 
vertical coordinate in the three dimensional mass transport 
equation such that there is no lateral variation of the 
total vertical length. This technique alJ.ows the use of an 
expanding and translating Eulerian grid system in the X-Y 
plane, moving within a time and spatially varying bottom 
and free surface. 

The complex grid system of Figure 3 . 1 in cartesian 
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space cani be transformed into an equal volume cell system 
in space, of Figure 3.2 . The transformation of the 

vertical coordinate proceeds from cartesian Z-space to 
dimensionless space by setting the the value of to 0, 
and -1 at the free surface and bottom respectively, 
(Figure 3.1) vhere 


Tl(x,y,z,t) 

_ z - 5 (x,y,t) 
H (x, y, t) 

(3.7) 

H(x, y, t) 

= h(x,y) + 5 (x, y, t) 

(3.8) 

T1(x, y, z. 

_ z - y, t) 

h(x,y) + § (x,y,t) 

(3.9) 


z 


= vertical coordinate (positive down) 


h(x,y) = depth at M.S.L. 


^(x, y, t) = free surface height above M.S.L. 

due to tidal variations 


The procedure is to non-dimensionalize the three 
dimensional mass transport Eq. 3.4 in the Z-coordinate to 
dimensionless ir^-space using the following transformation 
identities after Gordon (36) . The derivation of these 
identities may be found in Appendix A. 


- 5 - ( ) = 

H Til 


(3.10a) 


(3.10b) 
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> "eT < > - < > 


(3.10c) 


j4' + ( , O.iod, 


Equation 3.4 was split into seven terms and 
transformed using the transform identities 3.10 a,b,c,d. 

The mathematics and assumptions utilized in the 
transformation are treated in detail in Appendix B. for the 
interested reader. The transformed pseudo total velocity 
form may be written as 


^ + ± [U cl + i- [V cl + ^ [(JU cl = 0 
6t 6x T •' 6y T •* '• T 


(3.11) 


Where 0^, , CO^ , may be defined as 


= V 
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(3.12) 


(3.13) 


(3.14) 
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“s -ff (“s - ill'll 

+ V [('ll + 1) II + T1 ||l}} (3.15) 

Equation 3.11, written in flux conservative form, may now 
handle the transport and deposition of Lagrangian marker 
particles in X^Y-i^ space. 

The Eulerian three dimensional grid consists of equal 
volume cells throughout space, bounded at the free surface 
and bottom. The system may now take advective fluid 

velocities from an external source such as a two 
dimensional vertically averaged tidal hydrodynamic model, 
or a three dimensional tidal hydrodynamic model that 
utilizes the same transformation procedure (36). The SAPIC 
grid system can expand or translate in either the X or X 
dimension but there can be no expansion or translation in 
the ir^-diaiension due to the nature of the transformation. 
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IV. WAPIC ALGOBITHM 

The solution to the pseudo velocity equations is 
accomplished by a Laqrangian marker particle in Eulerian 
cell technique, utilizing Sklarcw’s total pseudo flux 
method. WAPIC employs a space staggered three dimensional 
Eulerian grid system composed of equal volume rectangular 
cells ( Figure 4.1 ) Concentrations are defined at the 

center of the cells and the advection velocities, U, V, w 
and the diffusion velocities, Uj-,, V]~, are defined at the 
center of the respective cell faces. Marker particles 
placed inside the Eulerian grid, which statistically 
represent the pollutant cloud are defined by their 
individual Lagrangian coordinates. A time cycle of the 
code can he divided into an Eulerian and Lagrangian step. 
Each step is explained as follows. 

aFPEODUGiBn^iTy- of the 

©BiOiNAL PAGE IS POOB 

A. Eulerian Step 

At the beginining of each time step, the 
concentrations, defined at ' the center of each cell, are 
used to calculate the diffusion Velocities defined at the 
center of the cell faces. WAPIC employs a centered space 
finite difference approximation for the diffusion 
velocities* The finite difference forms for the diffusion 
velocities may be written in cartesian space, or in the 
vertically transformed dimensionless space. 

1* Cartesian Finite Differences 
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The cartesian form of the partial differential 
equation of mass transport is used to model unbounded and 
uniformly bounded mass transport in three dimensional X-Y-2 
space. The partial differentials for the pseudo diffusion 
velocities ( Rgs. 1.6-a^b,c,) are approximated by centered 
space finite difference formulations 
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- c. 
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(4.1) 


(4.2) 


(4.3) 


The integer indicies, I, J, K, of the concentration field 

•r 

C, represent the three dimensional cell centers at which 
the concentration values are evaluated or defined* Since 
the concentration gradient in a given dimension is 
evaluated by the forward space approximation across cell 
centers from I to Itl, the diffusion velocity gradient is 
defined at the respective cell faces. These diffusion 
velocities are calculated throughout the the three 
dimensional space and then added to the externally supplied 
advection velocities, which also are defined at the grid 
faces, to yield the pseudo total velocities, Wy . 

2. Transformed Finite Differences 
As previously stated KAPtC atili>ZBS the dimensionlass 


I 


vertically transformed version of the pseudo total velocity 
equation, when simulating mass transport in a tidally 
doainated coastal -zone environment which has a time and 
spatially varying free surface and irregular bottom 
topography* The turbulent diffusion velocities in Bqs. 
3.12,13,14,15 are approximated by the forward space finite 
differences 
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As before, 
externally 
yield the 


these diffusion velocities 
supplied mass consistent 
total pseudo velocities. 


ate then added to the 
adyection field to 
Ot, Vt, Wt , in X-i>N 


space. 


B. Lagrangian step 

Following the calculation of the pseudo total 
velocities throughout the Eulerian grid, each marker 
particle is advected in three space by the velocity field 
from time T to T 4AT. Two methods of advance are employed 

1. Pure Explicit 

2. Time Cent'ered Iterative 

both using a bilinear weighting scheme to interpolate 
velocities for particles inside the grids. These 
procedures are as follows. 

1. Explicit Formulation 

The explicit formulation is used for simulation of 
mass transport of Lagrangian marker particles in steady 
undirectional or uniformly sheared flow fields (non curved 
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streamlines ) . Each marker particle within the Kulerian 
grid is transported from time T to T ♦ftT, a distance in 
X-Y-Z space UAt, VAt, WAt» respectively, respectively. The 
new Lagrangian coordinates become 


^EW “ ^OLD 


(4.8a) 


'NEW 


'OLD 


+ V_ 


At 


(4.8b) 


^NEW = ^OLD 

P P P 

where are the bilinear weighted velocities 

for each particle. The bilinear velocity weighting method 

is popularly used in the MAC method (16,27) for 

interpolating marker particle velocities in an Bulerian 

velocity grid. Figure 4.2 details a two dimensional 

p P P 

example for determining the , V^ , and velocities for a 
particle within the space staggered grid system. This 
procedure was generalised to obtain a three dimensional 
bilinear weighted average. A given WAPIC grid may be 
surrounded by a raarimura of 26 cells. ( Figure 4.3) The 
WAPIC solution determines the particular cell in which the 
particle's coordinates are located, then divides that cell 
into eight equal volume octants, and finds the octant in 
which the particle is located. { Figure 4,4) Each marker 
particle is assigned a "fictitious" volume equal to the 
volume of an Kulerian grid cell. This fictitious volume 
may overlap into a maximum of eight surrounding grid cells. 
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Thus a particle’s velocity is determined by a three 
dimensional weight of the nearest, eight, surrounding 
pseudo total velocities. Figure 4.5 details planar X-Y, 
X-Z, Y-Z views of the eight surrounding velocities used in 
the three dimensional bilinear scheme. 

2. Iterative Time Centered Approach 
The explicit particle movement technique with the 
bilinear weighting scheme was verified by Lange (28,29) 
against simple analytic diffusion cases in unbounded 
uniform and uniformly sheared flows to within 5.0 % 

accuracy of concentration prediction. When the streamlines 
of the flow field are curved, neutrally buoyant marker 
particles tend to develop trajectories not co-incident with 
streamlines of the flow field (37) . This interpolation 
error inherent in an explicit scheme with a bilinear 
velocity weighting technique, may be corrected by 
decreasing the calculational time step however at the 
expense of increased computational time. 

1 I 

Forester (37) proposed ah iterative time centered 
approach for prediction of marker particle motion in curved 
flow fields. Instead of employing equations 4.8a,b,c to 
predict particle movement, the following procedure is 
employed: 


NEW 


= X 


OLD 


+ (U, 


pOLD 


T 


+ U, 


pNEW 


) At/2 


( 4 . 9 ) 


where , represent the bilinear weighted 

particle velocities at the new positioTi. The method 







employs an iteration procedure where each distance estimate 
is used to produce* a bilinear weighted particle velocity at 
the new position. Using Eg. 4.9 and the new particle 

velocity at that position, a new position is calculated and 
is tested for convergence against the previous position. 

As an illustration of the utility of this technique, a 
two dimensional uniform circular vortex field was specified 
with a shear rate of du/dy= -0.0004 sec^ , and dv/dx= 
0.0004 seel with minimum and maximum velocities of -1.0 and 
1.0 ft/sec in each dimension respectively. The HAPIC grid 
consisted of a 10x10 X-Y matrix, thus specifying a time 
step of 250 for grid lengths in X and Y of 500 ft., using a 
Courant dynamic stability time step restriction for 
accuracy of particle motion prediction. ( Figure 4.6 ) 
When employing the non-iterative explicit bilinear weighted 
movement technique a massless marker particle placed in the 
circular vortex field was found to eventually spiral out of 
the field. ( Figure 4,7) When using Forester's time 
centered iterative approach, a particle placed in the same 
initial position in the vortex field, showed a precise 
circular trajectory. 

When the time centered iterative procedure is not 
employed for predicting particle transport in flow 
situatins with an eddy or vortex structure and a Courant 
time step limitation is used with the bilinear velocity 
interpolation scheme for particle motion, the deduction can 
be made that the mass or particles within the eddy will be 
accelerated out of the flow field prematurely. For 
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modelinq mass transport in complicated non uniform flow 
situations, the iterative time centered method requires 

: i 

much less time than a purely explicit technique for 
accurate prediction of particle motion. Exactly how much 
time will be saved for accuracy of particle motion 
prediction by using the iterative technique as opposed to 
the explicit' technique is highly problem dependent. 
Results by Spaulding (32) for accuracy of prediction of 
particle trajectory in a uniform circular vortex field 
field indicate a factor of 3 to 5 decrease in computational 
time costs with Forester's implicit technique as compared 
to a purely explicit particle movement technique with the 
same Courant accuracy • number . 


C. Expanding and Translating Grid 
The Lagrangian step contains the options of allowing 
selective expansion of the Eulerian grid system with the 
diffusive growth of the marker particle cloud and selective 
translation of the grid With the mean motion of the 
advected cloud, for any of the three dimensions in 
cartesian X-Y-Z space, and in the X-Y directions for the 
X-Y-ll vertically non-dimensional coordinate system. 
Incorporation of this procedure permits maximum resolution 
of the marker particle concentration field with a minimum 
number of computational grids. The expansion and 
translation procedures whi6h are particularly useful for 
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simulating instantaneous point releases in horizontally 
unbounded flow situations are developed in WAPIC as 
follows: 

The expanding grid system can operate in any 
combination of the three dimensions. During the expansion 
step each particle is tested for placement within a certain 
distance of an outer grid system boundary. Should a 

particle reach a grid cell boundary in any given dimension, 
then each grid cell expands a specified percentage of grid 
cell length in that given dimension. Figure 4.8 

illustrates an expansion in the X-direction for a 2x2 grid. 

As shown, the lower left hand corner of the grid is moved 
in the negative X-direction a distance one half of the grid 
cell length increase. WAPIC monitors the coordinates of 
this grid "origin" from time T=0, and updates its location 
in 3-space whenever there is an expansion or translation of 
the grid in succeeding time steps. 

The grid translation procedure is designed to follow 
the mean motion of the marher particle cloud as it advects 
through space. At time T, an average particle velocity 

(U , V , W ) is calculated in each dimension, from the 
P P P 

average values of the bilinear weighted particle 
velocities. The translation of the grid, in a given 
dimension, is achieved by moving each particle a distance U * 

ir 

^T from the old particle position. 
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and subsequently the origin { for example- in the X 
dimension ) is moved 


XR 


new 


= XR 


old 


+ u Zit 
p 


( 4 . 11 ) 


D. Boundary Conditions 

When the Eulerian grid no longer translates or expands 
in a given dimension, then boundary conditions must be 
employed in HAPIC for that given dimension. Because WAPIC 
is part Eulerian and part Lagrangian, boundary conditions 
are imposed on the Eulerian total pseudo velocity field, 
and on the Lagrangian marker particles within the field. 
Each three dimensional ceil is assigned a boundary 
identification number for logical control of the pseudo 
velocity calculations and movement of Lagrangian marker 
particles within the given cell. Table 4.1 summarizes the 
cell types or classif ications and the Eulerian and 
lagrangian boundary conditions imposed within the cells. 

1 . Eulerian Boundary Conditions 
As specified in Table 4 .1 , the Eulerian grid cells are 
classified into two groupings; these are - 

(a) Vertical boundary cells (reflection, 
deposition* and erosion) whose cell faces lie on a 
horizontal plane boundary 

(b) Horizontal cells (land, water, and open 
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CELL 

DESCRIPTION 

EULERIAN 

boundary CONDITION 

LAGRANGIAN 

boundary 

CONDITION 

LAND CELL 

NO FLUX 
^D == ° 

PARTICLE 

REFLECTION 

WATER CELL 

NONE 

NONE 

OPEN BOUNDARY 

o 

II 

PARTICLE 

REMOVAL 

VERTICAL REFLEC- 
TION 

^ or = 0 

D ' D 

PARTICLE 

REFLECTION 

VERTICAL DEPOSI- 
TION 

W ^ or W =0 
D D 

PARTICLE 

DEPOSITION 

VERTICAL EROSION 

SPECIFY BTM. SHEAR 
STRESS FOR EROSION 

PARTICLE 

ENTRAINMENT 


TABLE 4.1 Classification of Bulerian and Lagrangian Boundary Conditions 
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vertical plane boundaries) which are boundary and 
non-boundary cells not on a horizontal plane boundary. 

Through the use of the vertical transformation, WAPIC 
can treat the complicated bottom topography and time 
varying free surface of estuarine and coastal hydrodynamic 
flows as uniform boundaries within the space staggered 
finite difference grid. Rulerian boundary conditions are 
imposed by controlling the calculation of the diffusion 
velocities throughout the grid. Consider the construction 
of a three dimensional (fixed ) Eulerian grid with 
non-uniform lateral ( Y ) boundaries and a uniform free 
surtace and bottom, looking at a X-Z plane view. Figure 
4.9, it is seen that the upper row of vertical velocities 
directly coincide with the boundary. Thus their diffusion 
velocities would be identically zero by the definition of 
no flux across a land or solid boundary. The lower row of 
vertical velocities which do not lie directly on the 
boundary have finite diffusion velocities by the forward 
space finite difference approximation. Following Lange's 
(28) assumption of constant marker particle flux at an open 
boundary ( which rests on the assumption that the marker 
field is fully developed in the direction of the flow when 
it reaches the boundary ), the diffusion velocities are set 
identically equal to zero for the left hand column of 
u -velocit les . Figure 4. TO illustrates a horizontal plane 

view of the u-v diffusion velocity. For example, the 
diffusion velocities at 025 and 022 by the forward space 
finite difference Fq. 4.1 are zero. They do not have to be 
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equal to zero because the forward space concentration 
gradients can be finite by this definition, but are set to 
zero because they are located on the boundary. These 
assumptions are consistent with the no- flux assumption 
across a solid boundary. 

Each marker particle is treated as though possessing a 
fictitious volume equal to a grid cell, centered at its 
coordinates. Should this fictitious volume overlap onto a 
horizontal land or open boundary cell, then the particle is 
either reflected or removed from the grid field, ( Figure 
4.11 ) A marker particle may penetrate halfway into a 
horizontal land or open boundary before the horizontal 
boundary conditions applies. For the purpose of vertical 
resolution of marker particle movement, a particle may 
penetrate as far as the cell face on which the vertical 
cell boundary is defined. Particles reaching a depositonal 
vertical boundary are removed from the field by setting 
their Z-coordinates to zero and storing the X-Y deposition 
coordinates . 


E, Concentration Calculation 
The final segment of the Lagrangian step is the 
calculation ot the new concentration field from the new 
particle positions. Each marker is assigned a physical 
mass and a fictitious volUiBe the size of one Eulerian grid 
cell. For a given marker particle, the WAPIC code locates 
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the individual Eulerian cell and one of eight possible 
octants, of that cell, in which the particle's fictitious 
volume is located. Figure 4.12 details the concentration 
weighting scheme for a marker particle in a two dimensional 
grid system, WAPIC uses an extended version to allow for a 
third dimension. This mass sharing procedure is simply 
repeated for each particle in the field for calculation of 
the total cloud distribution of concentration in the grid 
system. 

F. Two Dimensional Velocity Processor 
The coupling of WAPIC to a two or three dimensional 
finite difference mesh Bulerian hydrodynamic velocity field 
requires the use of a "velocity processor". Since the 
WAPIC grid system is translating and expanding with the 
marker particle cloud as it advects with the tidal flow in 
a two dimensionally vertically averaged velocity field, the 
WAPIC Eulerian grid mesh points for velocities will not 
necessarily coincide with those of the externally supplied 
fixed Eulerian velocity field. Following the work of 
Dobson (38) , a least squares surface fitting technique has 
been adopted to interpolate u and V horizontal velocities 
from a two dimensional space staggered finite difference 
grid system. The method employs the twelve nearest fixed 
grid neighbors (in the form of a cross ) around the point 
of interest , and fits a quadratic surface of the type - 

2 2 

~ ®2 ^ ®3 ^ ( 4 . 12 ) 

e = Matrix Coefficients 





4* 

|-* — M- 


+Cl+l,j+l 

4 - 


Cl+l,J 


v-CELL VOLUME 
(AXAYAZ) 

Mi-PARTICLOMASS 

C-CONCENTRATION 


PARTICLE 


CONCENTRATION COMPUTATIONS 

Cj^ j = ^ UZ) - XP) - YP) 


^I^ l .J = ^ (AZ) (XP - - VP) 

= 72 (^mxp - ^ - ^ 

- XPKVP - 


FIGURE 4.12 Dimensional X-y View of rarticlc Goncon- 

tration Weighting Scheme. Inside tl'jc 
Eulcrian Grid Matrix ... .... .... 


51 . 


The coefficients of this eguation are deter mined by an 
analytic matrix solution procedure. ( Figure 4.13) The 
velocity processor interpolates velocities from the fixed 
Eulerian system to the WAPIC mesh at each time step of the 
tidally driven transport simulation. This velocity 
int erpcla tion technigue is ideal when the WAPIC total grid 
mesh span is smaller than two or three hydrodynamic grid 
cells. When the WAPIC Eulerian grid grows to such a size 
that it is as large as a 4x4 tidal hydrodynamic field, then 
the surface fitting technique may produce discontinuities 
in the interpolated velocity field, if there are large 
variations in the velocity vector field. Should the 
expanding grid system encounter an island or outer 
boundary, expansion and translation is stopped and the 
twelve point surface fitting technique is replaced by a 
simple linear velocity extrapolation technique using the 
four nearest hydrodynaraic grid neighbors to interpolate the 
velocity at the WAPIC grid point. 

G. Particle Generation Procedures 
WAPIC has three methods for locating laarker particles 
within the Eulerian cells. These methods are as follows: 

1, instantaneous Gaussian Belease 
In general, turbulent *Ficlcian diffusion from a point 
source in homogeneous uni'-directional unbounded flow has a 
three dimensional Gaussian character (33). The 
concentration field can be described as exhibiting normal 
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behavior around a maxirauni concentration. The assumption is 
made in the WAPIC procedure that the so called 
instantaneous point release, after a sufficient dispersion 
time, has some Gaussian character which can be described by 
its standard deviations about some mean position in 
cartesian space. The assumption is made that an initially 
finitely sized evenly distributed waste cloud such as the 
turbid plume from a dredge spoil disposal operaton can be 
represented in the long term dispersion phase as if it were 
a point release that had grown to a normal ” Gaussian" 
concentration distribution in space. 

Several experiments were conducted using Gaussian 
random number generators* similar to those used by Lange 
(28) , to place the marker particles inside the Rulerian 
grid cells. These techniques proved to give particle 
distributions that displayed visual "bunching" or crowding 
around the major area of release center and required 
excessive computational time to place a large number of 
marker particle within the grid system. This problem was 
directly related to the quality of the Gaussian random 
number generators. 

An alternative method was devised. This method 
specifies the number of particles to be placed in each ceil 
by using a Gaussian probability function. The assumption 
is made that 99.5% of the particles are within a range of 
i 3+sigma ( sigma= std. deviati-on from a given mean in 
three dimensional space ). By specifying the three 
dimensional grid cell and standard deviation lengths and 
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the distance of each Eulerian grid cell from the mean 

particle position, the total number of particles per 

i r 

Eulerian <jrid cell can be determined by a simple analytic 
Gaussian function. A uniform random number generator is 
used (3.9) to locate the appropriate number of particles in 
each Eulerian grid cell. Should the need arise to utilize 
another type of initial distribution of marker particles, 
all that is necessary is to specify another analytic 
function that will specify the number of particles per 
Bulerian cell, 

2. Continuous Gaussian Releases 

A plume or continuous Gaussian release may be 
simulated by releasing a Gaussian particle distribution at 
the point source location for each succesive time step. 

3. Uniform Random Distribution 

The option has been included in WAPIC to fill any 
rectangular volume in the Eulerian grid with a uniform 
random marker particle distribution. This method may be 
useful in simulating the short term transport from an 
instantaneous plug or plane constituent release. 


H. Particle Settling Velocities 
The advantage of utilizing the PIC method for modeling 
transport of suspended sediment particles is that a 
spectrum of settling rates for particles within the cloud 
may be superimposed upon the advective and diffusive 
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particle velocities. The assamptions are made, for sake of 
simplicity and due to the lack of adequate field 
measurements of vertical density, that the water column is 
vertically well mixed in density and homogeneous with 
respect to vertical turbulence; and that a given Lagrangian 
marker particle statistically represents a given mass of 
sediment particles in the water column, all of equivalent 
density and still water asymptotic settling rate, Hs. 

This particle settling rate is calculated by Watson's 
(40) empirically modified version of Wubey's (41) 
analytically derived still water settling law for sand 
grains. 



+ (4/3)SR^(p - 

P 


SRp 


FL 


g - 3Ap.) 



(4.13) 


S = pressure drag coefficient = 0.5303 dimensionless 

A = viscous drag coefficient = 0.623 dimensionless 

R = particle radius cm 

W = settling velocity cm/sec 

p - fluid dynamic viscosity poise 

3 

P - mass density gm/cm 

PL, P = subscripts for fluid or particle density 
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This vertical settling velocity is then added to the 
vertical advective and diffusive velocities to yield 

W = W + (4.14) 

T D S 


or in dimensionless vertical velocity form 

(D = (JU + (4.15) 

T D S 


w_ W 

T 

D S 



Using 

the 

above equation. 

a spectrum of particle 

settling 

rates 

can 

be represented 

by specifying sediment 

particle 


diameters and densities for the individual Lagrangian 
marker particles. 


57 . 


V. ACCUARCy PARAMETERS, PARTICLE-CELL RELATIONSHIPS 
Since WAPIC is a Lagrangian particle in Eulerian 
finite difference cell technique, the questions may be 
asked: What effects do the relationships between the 

number of marker particles and the number of Eulerian cells 
used to resolve their distribution have on the accuracy and 
stability of WAPIC'S predictions? How well does HAPIC'S 
three dimensional Gaussian distribution procedure produce a 
truly normal cloud distribution? A discussion of the 
particle generation procedures and the accuracy and 
stability requirements of the particle in cell method 
follows. 


A. Number of Eulerian Cells 

A basic requirement for all methods that use Eulerian 
grid systems and finite difference or finite element 
procedures is that sufficient grids are utilized to resolve 
the concentration gradients of interest. In order to make 
an appropriate estimate of the number of Eulerian grid 
cells needed to resolve the normal concentration 
distribution in time and space, the technique employed by 
Lange (28) will be used. 

Consider the one dimensional finite difference 


a pptoxima tion to the diffusion velocity 

' k: 

u 




(5,1) 


i + + 1/2 Ax 

by dropping the subscript for the i+1/2 terms and expanding 
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Ci+1 and Ci in a Taylor series about i-t-1/2 to the first 
contributing terra; 


series 


A X + i iLS. A 

^ix “ 24 ^ 3 ** 


X + H.0.T.3 


(5.2) 


Expanding Kg, 5.2 and dropping the higher order terms 


yields 


K ^ K . 2 .3 

JX #c X ( Xx J|_C 

"c” ^fx ”24 c ■ ^^3 


(5.3) 


The first terra on the right hand side of Rq. 5.3 is 


the exact differential expression for the diffusivity 
velocity, and the second terra is the first contributing 


error terra. Taxing the ratio of the WAPIC diffusion 
velocity approximation and the exact differential diffusion 
velocity yields 


g' .D. = T + / « AiL. 

“exkt >=■ 


(5.4) 


Eulerian grid cell size on the diffusion velocity error 
terra in Eg. 5.4, consider the simulation of a one 
diraensional Gaussian concentration distribution 


Or 2 ^2^ 
= ^ L -x /2a 3 


(5.5) 


CT = standard deviation 


0 = Mass released 
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The first and third derivatives of F.q. 5,4 can be defined 
as 




Qx r 2 2- 

^ e [-X /2a 1 


(5.6) 


Qx r 

5 ^ 7 


2 o o 

e [-x^/2a^^ 


(5.7) 


Substitution of Eqs. 5.6 and 5.7 into Eq, 5.4 , yields 


u. 


u. 


F.D. _ ^ 


EXACT 




(5.8) 


Since 9S% of all particles in a one dimensional Gaussian 
distribution lie within a distance x=3^sigKa ,the 
substitution of x=3*sigoa into Eg. 5.8, yields 
u. . 2 ■ - 

_ = 1 : (5.9) 

U 2 

EXACT 4a 

For purpose of illustration, the substitution of 
sigma- AX or 2* AX in Eg. 5, 9 yields a diffusion velocity 
error of 25?> and 6.25 1 respectively. This indicates that 
as more grid cells are used to resolve a given standard 
deviation of the Gaussian particle cloud, the WAPIC finite 
difference approximation becomes an increasingly better 
estimate of the exact analytic diffusion velocity. 
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B. Number o£ particles Per Cell 

The Laqrangian marker particles statistically 
represent the dissolvecl or suspended mass in the ambient 
flow field. Usually, for the case of numerical modeling of 
horizontal mass transport in estuarine and coastal areas, 
the turbulent diffusive transport is on the order of 0 . 1-5 
% of the advective fluid velocities ( 35 ). When the 

particle cloud is small in relation to the scale of spatial 
variations in the advective velocity field, then the 
diffusive transport will govern the spread of the cloud 
about its center of mass, while the center of mass being 
transported with the mean velocity of the flow field. When 
the particulate cloud reaches into a field of velocity 
shear, then the advective velocities will dominate the 
spread of the particle cloud. Therefore, the more 
particles that can be used to represent a cloud or mass 
distribution in space, then the more finely can the 
diffusive mass transport process be simulated, provided 
that there is accurate resolution of the velocity field. 
However computer storage ceguirements and run time ( run 
time is directly proportional to the number of marker 
particles ) may place limits on the number of particles 
available for a given simulation. 

WAPIC stores { IBB 370-155 , Fortran IV ) the X,y,Z 
coordinates and mass of each particle in one dimensioual 
arrays. Figure 5.1 illustrates the minimum amount of 
storage ( real*4 ) bytes required to handle these arrays 
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versus nujrber of particles. 

Assuming that the marker particle cloud is three 
dimensionally Gaussian, then the length or range of the 
cloud is approximately h^sigma ( 3^sigma around the mean, 
where sigma is the standard deviation of the distribution) . 
A ratio of sigma/dx of 2, ( where sigma/dx equals the 

standard deviation to grid length ratio ) provides a 
diffusion velocity with an approximate error of 6 . 

Thus, if 6 % error is acceptable, there will be required at 
least 12 grids (6x2) to represent the particle cloud in 
each dimension. Assume then that a 15x15x15 X-Y-Z Eulerian 
grid can resolve the concentration gradients within the 
Gaussian cloud, WAPIC utilizes three dimensional storage 
arrays for the Eulerian U, V, W velocities and 
concentration C and boundary cell flag M. Thus a 15x15x15 
grid requires at least another 

5x15x15x15x4 = 67.5 k bytes of real*4 

Using Figure 5.1 and the above, for a Gaussian cloud 
of 5000 marker particles, at least 157.5 k storage bytes 
are required to accurately predict the three dimensional 
diffusion velocity. 

As explained in Chapter IV. , WAPIC assigns a 
fictitious cell volume the size of an Eulerian grid to each 
marker particle. This technique allows the overlap of 
particle mass into its eight nearest neighboring cells, 
which tends to smooth the concentration field and thus the 
diffusion velocities. The question might be asked - How 
few marker particles per Eulerian cell, when using the 
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fictitious volume sharing scheme, will yield meaningful 
results for particle diffusion velocity calculations and 
therefore the diffusive spread of the particle cloud ? 
Lange (28) determined empirically that - ” As few as one 
particle will obtain meaningful results. With fewer 
particles than that, ie. when particles have no neighbors 
around them, the diffusion velocity algorithm moves the 
particle to a center of an Eulerian cell and freezes it 
there , 


C, Stability, Accuracy , Time Step Restrictions 
Since WAPIC uses a finite difference Eulerian grid, it 
has a dynamic stability or time step restriction that 
states that the fastest moving particle can translate no 
more than one half a cell length in a given time step dt. 
This restriction is commonly used in vorticity transport 
solution techniques (16) and numerical tidal hydrodynamic 
models (17,18,42), and is defined for WAPIC’S uses as 


^ AX ' aV 


Az 



1 

2 


( 5 . 10 ) 


where 0^ , , W^are the Eulerian total pseudo transport 

velocities at time T, and AX, Ay , Az are the grid cell 
lengths. For purposes of this work the above time step 
restriction in Eg. 5.10 shall be called an accuracy 
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requirement, or accuracy number. Roache (16) states that 
PIC methods exhibit high frequency oscillations in the 
marker particle field which may be visualized as a typo of 
instability in the solution procedure. The time step 
restriction acts to damp the amplitude of these marker 
particle oscillations in the Eulerian grid space. An 
example of this solution instability and the effect of the 
time step restriction will be demonstrated in a subsequent 
chapter. 

As discussed in Chapter IV. ^ should the user choose 
not to utilize the time centered iterative particle 
velocity interpolation technique and the flow field 
displays an eddy structure or curved streamlines, then the 
time step should be reduced below the value given by Eq. 
5.10, for accuracy of particle motion prediction. The 
exact value of the time step needed is situation specific. 
The author recommends, though, that the time centered 
iterative technique be utilized, because of the substantial 
savings in reduced computational time for a specified 
accuracy of particle motion prediction (37), 


D, Choosing a Standard Deviation Grid Size Ratio 
As discussed in Chapter IV,, WAPIC simulates an 
instantaneous Gaussian release in space by using a normal 
distribution equation to specify the number of particles 
per Eulerian cell. A uniform random number generator 

THE 
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locates the individual particles in each Eulerian cell. 
The question was asked how well does this procedure produce 
a truly normal particle distribution? This question arises 
from the knowledge that the initial particle distribution 
influences the short and intermediate time concentration 
field when modeling mass transport in homogeneous turbulent 
unidirectional flow (33) . Also, as mentioned previously, 
the ratio of the initial cloud standard deviation and grid 
spacing determines the accuracy of WAPIC’S diffusion 
velocity prediction. A rephrased question might be - How 
many marker particles are necessary to achieve a normal 
distribution with a given sigma/dx, ( where again siqraa/dx 
represents the ratio of the standard deviation to grid 
length ), ratio that will allow accurate short and long 
term prediction of cloud dispersion ? 

In order to achieve an accurate diffusion velocity and 
thus an intermediate and short terra cloud disperson 
prediction, a sigraa/dx ratio greater than one should be 
used, and preferably a ratio of two, which would be 94% 
accurate in the diffusion velocity approximation. As the 
ratio of sigma/dr increases beyond two, the computer 
storage and run time increase because the size of the 
Eulerian grid must increase. As in pure Eulerian schemes, 
a trade off exists between accuracy of prediction and 
computational costs, 

Yevgevich. (43) *• provides a method for testing a 

particle distribution for normality. He provides two test 
criteria for a Gaussian or normal character 
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-0.10 < Cs < 0. 10 


and 


-0.50 < E < 0.5 

where Cs is the sXewness coefficient and E is the excess 
coefficient of the particle distribution. The skewness 
coefficient is a measure of the asymmetry independent of 
the dimension of the variable, and is defined as the ratio 
of the third central moment to the second central moment to 
the 3/2 power. 


M. 


M, 


C = 
s 


372 

M2 


= Second Moment 


= Third Moment 


C = Skevmess Coefficient 
s 


(5.11) 


where 


„ 

r N 


N 

£ 

i=-l 




- X) 


X = Mean 
N = # Samples 


M = Moment # 
r 


r =1, 2, 3.. 


(5.12) 


The excess coefficient is defined as the excess kurtosis 
relative to that of the normal distribution 


E = Y 2 - 3 , . 

Y2 = ^4/^2^ = ^4/^7^ 


(5.13) 


(5.14) 






where a normal distribution is defined as having a kurtosis 
coefficient of 3. If E is negative then the distribution 
has a sharper peak than the normal distribution, and if E 
is positive then the inverse is true. 

A procedure was written to test MAPIC'S Gaussian 
particle distribution scheme for ability to produce a 
"truly*' normal distribution for varying; numbers of 
particles at various signia/dx ratios. Three dimensional 
"spherically normal" particle distributions were produced 
for siqma/dx ratios of 1.0 to 3.0 iu increments of 0.25, 
using 1000 to 10000 particles in increments of 1000 for 
each ratio. The second, third and fourth moments around 
the mean were calculated, which allowed the calculation of 
the actual sigma/dx ratio, skewness coefficient and excess 
coefficient of each distribution. 

Figure 5.2 illustrates the actual sigma/dx ratios 
derived from the estimated sigma/dx ratios at various 
particle numbers. From this figure the following 
approximations can be made. 

1. Approximately 1000,1700,3000,4500 particles 
are needed to achieve an actual sigma/dx ratio of 
1.0,1.25,1.5 and 1.75 respectively. 

2. At estimated ratios greater than or equal to 
2.0 more than 10000 particles are necessary to achieve an 
exact ratio agreement with the estimate* 

3. At an estimated ratio of 2.0, there is little 
incentive to use more than 3000 particles, because there is 
very little change of slope in the actual line from 3000 to 



I' 



RCTURL SIGhRX/DX 


MRRKER PRRTICLE SIMULRTION OF 3-D GRUSSIRN DISTRIBUTION 
RRTIO CLOUD X-STRNDRRD DEVIRTION & GRID SPRCING VS PRRTICLES 


FOR ESTIMRTES OF SIGMRX/DX RRTIOS 



10.00 TO. 00 30. 00 40. 00 50. 00 60. 00 70. 00 80. 00 30. 00 100. 00 


NUMBER OF PRRTICLES lo^t 2 

FIGURE 5.2 Actual Standard Deviation of WAPIC Particle Distribution Versus Number of Particles for 
Estimated standard Deviation to Grid Length Ratio 




cr> 

CO 



69 . 


10000 particles. 

4, In order to achieve sigma/dx ratios greater 
than 2.0 within an accuracy of 10% more than 2000 particles 
are required . 

Figure 5.3 illustrates a plot of the particle 
distribution X-skewness coefficients at sigma/dx ratios of 
1.0, 2. 0,3.0 versus the number of marker particles. From 
this figure the conclusion is made that WAPIC'S procedure 
produces particle distributions well within the skewness 
constraints regardless of the number of particles and 
sigma/dx ratio, within the limits tested in this example. 

Figure 5.4 illustrates a plot of the X-excess 
coefficients at sigma/dx ratios of 1,0, 2. 0,3.0 versus the 
number of marker particles. Prom this figure the following 
conclusions are drawn : 

1. At sigma/dx ratio of 3.0, a negative excess 
kurtosis is produced when less than 3000 particles are 
used . 

2. WAPIC'S procedure produces an excess 
coefficient for a normal distribution for sigma/dx ratios 
of 1.0 and 2.0, well within the range of excess of a normal 
distribution for sigma/dx ratios of 1. 0 and 2.0 for all 
particle cases tested. 

The method produces a normal distribution with a 
slight negative kurtosis (greater than -0.5 ) for particie 
numbers greater than 3000, indicating slightly higher 
concentration gradients in the center of the distribution. 

In general. Figs. 5.3 and 5.4 illustrate the 
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effectiveness of WAPIC'S Gaussian particle distribution 
procedure in producing a "truly normal" distribution. 
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VI. TURBriLENT DIFFUSION AND DISPERSION 
A. Definition of Diffusion and Dispersion 

Within the literature of numerical modeling of mass 
transport the processes responsible for net transport of 
dissolved or suspended constituents in a turbulent fluid, 
are advectionj diffusion, and dispersion (4). Advection 
refers to the. transport of some dissolved or suspended 
constituent by a spatially or time averaged velocity 
component. 

A qualitative definition of diffusion and dispersion 
has been suggested by Eiolley (47) - ” Let diffusion refer 
to transport in a given direction due to the difference 
between the true convection in that direction and the time 
average in that direction. Let dispersion refer to the 

transport in a given direction due to the uncertainty in a 
given direction due to the spatial average in that 
direction,” 

The turbulent ensemble time averaging of the mass 
transport equation produces a mean advective flux term and 
non linear advective terms of the form (4,17) 


(U^' G') 


i = 1, 2, 3 


(6,1) 


This nonlinear advective term is impossible to specify at 
every given time and spatial point in a continuum. The 
assumption, which has physical validity, is that turbulent 
diffusion develops a ’• Gaussian” or normal concentration 


74 . 


distribution in space centered aroiincl a 
given dimension of measurement (33). 
assumption is made that 


mean value in that 
The linearizing 


1 _ = J_ K 1^- 

Jx, (Oj^' C) '"i Jx. 


where 


(U^' C) 


( 6 . 2 ) 


( 6 . 3 ) 


This allows the formidable time averaged values of the 
product of the deviations from the mean of the 
concentration and velocity field, to be replaced by the 
product of a constant coefficient and the gradient of the 
mean concentration. 

In a similar manner, the three dimensional mass 
transport equation may be also spatially averaged with 
reference to deviations from the spatial mean velocity and 
concentrations. The gradients of the spatial average of 
the product of the deviations of concentration and 
velocities from the mean values are replaced by the product 
of an asymptotic ''dispersion” coefficient and the spatially 
averaged mean concentration gradient. 

The processes of turbulent diffusion and dispersion 
can be visualized as the "advective” transport by 
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compoiKints of velocity which deviate fccm the time aad 
spatially averaged mean component of the velocity, ie. 
'•tiirbulent velocity fluctuations.” Assuming that the 
molecular diffusion scale is small compared to the scale of 
"turbulent” motions of a fluid particle, turbulent 
diffusion and dispersion do not exist; there is only 
advection of suspended mass by random velocity components 
(33). The terms turbulent diffusion and dispersion are 
mathematical descriptions for extra terms left over from 
the time and spatial averaging of the three dimensional 
mass transport equation. 

Deterministic numerical mass transport models linearly 
simulate the turbulent transport of a suspended constituent 
through a distribution of finite grid cells or elements in 
space. Due to computer storage and computational time 
limitations, the numerical solution to the partial 
differential equation must be represented by discrete time 
steps across finite grid lengths. The continuum technique 
is to advance the initial concentration distribution 
through time across the finitely divided space domain. 
Numerical models thus employ both time and spatial 
averaging techniques to estimate a given level of 
con stituent . 

For the numerical modeling purposes of this work, 
Fickian diffusion and dispersion will be considered to be 
one and the same, as opposed to the definitions of Hoirey. 
in lieu of a more proper choice, the term turbulent 
^ diffusion will be used to describe the uncertainities in 
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the spatial and time averaged values of the mean advective 
and concentration field. 

Thus K is a coefficient relating the uncertainty in 
the measurement of the advective velocity and concentration 
field due to time or spatial averaging the instantaneous 
flow and concentration field. Components of motion which 
have length and time scales smaller than that used by the 
numerical or analytic solution technique, may appear as as 
sub grid turbulence or may not appear at all. Turbulent 
advective motion with time and length scales greater than 
that used in the numerical system are difficult to specify 
within the field, and as will be discussed shortly, are 
approximated empirically by a scale dependent diffusion 
law » 

Before presenting the methods tor estimating diffusion 
coefficients , a discussion of ” FicXian** diffusion and the 
"scale effect" in turbulent mass transfer would be 
informative. 


B. Fichian Diffusion 

As previously stated, the turbulent diffusion 
coefficient is used to replace the non linear terms of the 
form d/dx{ ), which result from the averaging process. 

The popularly termed Reynold’s shear stress analogy, 
employed by Yon Karraann in 1934 (4B) for a stationary 

homogeneous turbulent flow field, states that the rate of 
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transport of any characteristic of a fluid across a 
streamline ( heat, salinity, silt, etc. ) is proportional 
to the gradient of moaentura { or of the given property ) in 
the direction of transfer. 


u. ' u. ' 

TuJTFT 

This analogy allows turbulent diffusion to be defined in 
the same manner as molecular " Fickian” diffusion which has 
a *' Gaussian*' character (3 3). 


C. Statistical Properties of Diffusion 
The asymptotic Fickian description of turbulent eddy 
diffusivity assumes that the free stream turbulence has a 
stationary and homogeneous character. In the ocean, 
turbulence is not homogeneous or stationary. The 
observation was made (33), that eddies larger than the 
particulate cloud, advect it as whole, while smaller scale 
eddies produce a scattering of the cloud around its center 
of mass. As the cloud expands in size the boundaries 
between the eddies contributing to advection and those 
responsible for diffusion shift to larger scales. The 
cloud then goes through a short explosive expansion stage 
(33). The turbulent Fickian approach for describing cloud 


U. ' c' 

_i__ = K. = E. 

XI 


(6.4) 
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expansion as a linear function of time fails to explain 
this behavior. 

An alternate approach for describing the diffusion 
mechanism is the "stochastic" or statistical approach 
employed by Taylor (49) . Taylor's equation for one 
dimensional diffusion in a turbulent stationary homogeneous 
mean flow field is defined as 




■'u, . , , 
1 (t) 


(6.5) 


where ■ ^i(t) - the radius of inertia of the particle 

cloud, equal to the rms of instantaneous particle 
displacement ( x^) ^ and, 

R„» = u. * (t)' U. + T)/U;'^ (6.6) 

D • . . .X. .i. 

r.(T) 

is the Laqrangian correlation coefficient of the particle 
velocity distributions taken over the interval T to T + AT; 
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where 

- the time lag 

I ^ 

U* = the mean squared instantaneous turbulent velocity in 
the i-direction. 

Equation 6.5 describes the diffusion of a fluid 
particle about the mean position of the group. If the free 
stream turbulence is convected at a velocity Ui, then 
diffusion refers to the instantaneous displacement from a 
point moving with velocity IJi. a A general form of 61 ft) 
for all time t is difficult if not impossible obtain 
because of the functional dependence on the Lagrangian 
correlation coefficient. The Lagrangian correlation 
coefficient of the velocity field is inde pendent. of time 
for a stationary process and is an even function. This 
variable describes the "persistence’* of the value once it 
is realized. 

Equation 6.5 can be approximated, though, at the 
limits, T=0 and T-^infinity. At a very small diffusion 
time'tf-^O, where the lira Rnl (*t0 * T r 6,5 becomes 

(t) - u. (6.7) 

This relationship states that at small diffusion times, the 
standard deviation of the particle cloud grows 
approximately linearly with time. This phase is Called the 
linear phase, where the center of mass of the relatively 
small Cloud meanders, with reference to a fixed frame, with 
the turbulent motion of the eddies. The cloud is so small 
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that the relative diffusion is due to sub-inertial range 
eddies, 

When ^ infinity, the so called asymptotic phase of 
diffusion, the cloud is large compared to the energy 
containing eddies. The center of gravity is no longer 
subject to displacement by turbulent eddies and is 
translated solely by the mean motion, the turbulent 
diffusion coefficient being defined as 

h t = K. . (6.8) 

2 dt Lx Asymptotxc 

= Lagrangian time scale 



where K; , is a ” Ficlcian" constant diffusion coefficient, 
and the standard deviation of the cloud grows as the square 
root of time. 

During the intermediate stage of diffusion, 0 <^ < 
OO, the relative diffusion is dominated by eddies on the 
order of the size of the cloud. Dimensional considerations 
by Batchelor (50,51) , show that 


1 _ 1/3 4/3 

2 dt ® 


= K. (t) 

X 


Intermediate 


(6.9) 


Equation 6.9 states that the apparent eddy diffusivity is a 
non-linear function of time and increases with the 4/3 
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power of the cloud standard deviation. The validity of 
this relationship was first demonstrated by Stommel (52) . 
Eg. 6,9 may be integrated to yield (33) 

^ 2/3 ^ gl/3 ^ j (6.10) 

1 o 

In terms of time elapsed from the effective origin, t , the 
particulate cloud grows, or explodes, as (t - t^ )^* 

This behavior, is in contrast, to the (t)*'^*’ growth in the 
asymptotic Fickian diffusion stage, or in the initial 
linear growth stage. Csanady shows that the lifetime of 
this explosive stage is relatively short, or shorter than 
the lifetime of the energy containing eddies. As the cloud 
continues to grow the asymptotic diffusion stage is 
reached, where the cloud becomes large compared to the 
turbulent eddies. Csanady summarizes the governing 
mechanisms for horizontal cloud spread in the ocean: 

•’ If the release of marked fluid takes place some distance 
from the shore, the field of turbulence surrounding the 
effluent may be expected to be horizontally homogeneous to 
a good approximation. Our theoretical results should 

therefore apply to such releases; in particular, if the 
initial cloud size is small compared to the typical 
horizontal eddy, a phase of accelerated relative (4/3 law) 
diffusion should be observable, during which the cloud 
grows non-linearly with time. The duration of this phase 
should be relatively short, and be followed by the 
asymptotic phase, where the cloud grows with the square 


82 . 


root of time, •' 

Some striking support can be found for Csanady's 
observations in Sayre and Chang’s work (53). They cite the 
experimental evidence of Yotasukura, et. al. (5h) , Glover 
(55), and Godfrey and Frederick (56), for longitudinal 
diffusion in homogeneous turbulent bounded channel flow: 

" In general the agreement for Fickian theory and 
experimental observations is poor in the early stages of 
diffusion, but tends to improve with increasing diffusion 
time or distance from the source.” No comment was made 
about any attempt by Sayre and Chang to fit the initial 
data with a 4/3 type cloud growth law, though. They 
recommend that ” The Gaussian or asymptotic turbulent 
diffusion mechanism be used for predicting diffusion at 
large tines or distances from the source in flume diffusion 
studies. It 

What these flume diffusion experiments could not 
ascertain is that when a dye patch becomes on the order of 
1 kilometer or larger its diffusion rate increases sharply 
again, as if the 4/3 power was valid over a much larger 
range of eddy motion. Hunter and Spaulding (7) provido a 
review of attempts to validate the 4/3 law in the ocean, 
comparing the work of Yudelson (57) , Okubo (5fl) and Koh and 
Chang (5). Figure 6.1 shows a comparison of Okuho’s data 
and a least sguares fit (composite) for diffusion 
coefficient versus iength scale. Where the 4/3 law is 

K. = e ,1 i = X or y (6.11) 

1 IX 
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Hunter offers 

0 9983 ' 

= 0.10 1^ ' in cm. (6.12a) 


0 

= 0.0059 1^ * in feet (6.12b) 

This indicates a nearly a linear growth dependence for the 
diffusion coefficient with the Cloud length scale. Indeed, 
a theoretical point source release which expands to the 
order of kilometers would experience a cascade of turbulent 
eddy sizes or scales of motion thus experiencing a 
multitude of linear asymptotic, and non-linear scale 
dependent diffusion phases. The practical utility of 

fitting a straight line to the above data obscures the 
cascade of cloud growth phases. 

The discrepancies between simple statistical theory 
and actual diffusion theory are evident from the fact that 
turbulence in a flow field is not homogeneous, ie., many 
scales of eddy motion are present. As a particulate cloud 
expands, large turbulent eddy scales come into play. A 
constant eddy diffusivity plateau (asymptotic) is 
int erai ttently reached, and then surpassed as the cloud 
reaches into a larger scale flow structure, the cloud 
seeming to explode, growing with (t-t^. ) , for a 

relatively short time. 
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If these larger scale flow structures can be regarded 
as a velocity shear of the mean flow field, then the 
increased cloud growth rates may be attributed to the 
non-uniformity of this field (33). As a practical measure 
Csanady recommends that the homogeneous turbulent mean 
velocity field model be retained and the use of empirical 
( h/3 law) data on cloud growth be used to predict 

diffusion coefficients. 


D. Practical Calculation Methods 

For practical use, in predicting transport of 
constituents in numerical models, diffusion might now be 
defined as the turbulent deviations in the mean velocity 
field due to the time and spatial averaging processes. 

1. Horizontal Diffusion 

The wort of Elder (59) has been applied by Holley, 
Harleraan, and Fischer (60) to estimate horizontal diffusion 
coefficients in the tidally dominated estuarine 
environment. Hunter and Spaulding adopted F,lder*s method 
with the addition of Leendertse's (61) formulation for the 
added effects on diffusion due to the wind. 

K. =K. (Elder) (6.13a) 

where the Elder diffusion coefficient relations approximate 
the observed spatial variations in the mean flow field in a 
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given dimension (shear) to the turbulent diffusion 
coefficients ( popularly termed dispersion relations) 

K - 5.93HU* “ K = K Offshore (6.13b) 

X X y 

K = 0.23HU* = K Estuary (6.13c) 

y y 

K = 0.06HU* (6.13d) 

z 

and the added effects on the total diffusive transport due 
to the wind are approximated by a scale dependent 4/1 cloud 
growth law 

0 9983 

1 = 3 V^GT c ^ = 1 (6.13f) 

X 1 X y y 

These modifications approximate the added effects of the 
wind in the form of turbulent eddies, surface drift 

currents, and predicts diffusion when the mean velocity in 
the tidally driven flow field is zero. A first order 
estimate of the length scale of the diffusing cloud might 
be the physical horizontal Fulerian grid size, or a typical 
vertical depth of the hydrodYnamic field. Thus for a grid 
size equal to a nautical mile, ( 6800 ft. ), a first guess 
for an asymptotic diffusion coefficient might be 

0 9987 2 

K - 0.0059 (6076) * =40 Ft /sec (6.14) 

The horizontal diffusion terms can then be calculated at 
points throughout the Fulenain grid using Flder's 
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relationships, ( Eqs. 1 1 b,c,d ) , modified by the addition 
of a 4/3 power law to take in the added effects of the 
wind, { Eq . 6.13 e ), as large scale turbulence within the 
flow field. 

2. Vertical Diffusion and Sediment Transport 
As stated by Hunter and Spaulding (7) - " Diffusion in 
the vertical direction is restrained by the bottom and the 
free surface, and by the magnitude of the vertical 
component of the flow, which is normally very small 
relative to the horizontal flow. Thus it has an effect 
several orders of magnitude smaller than that of horizontal 
diffusion. The latter, however, is still a small effect 
when diffusive transport is compared with advective 
transport. This is not likely to be the case in the 
vertical direction, since the vertical flow is also very 
small. Vertical diffusion and vertical advection are very 

closely related effects, since they are driven by the same 
force (instability) , Since vertical diffusion is also 
highly dependent on wave action, it is likely to be more 
important than advection." Numerical mass transport 
schemes normally use diffusion alone to model vertical 
transport of constituents, because of ignorance of the 
vertical advective field. 

The procedure is to calculate the vertical diffusion 
coefficient by a correlation with the Richardson number (7) 

R. g'/p where p = density of water 

/ z = vertical level 

h = total depth 
U = horizontal velocity 


(6.15) 
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which relates the stabilizing forces of the vertical 
stratification and the destabilizing forces of the lateral 
shear flow. Table 6. 1 details a summary of the 
experimental values of the vertical diffusion coefficient 
in the ocean and Table 6.2 (from Spaulding (32)) details 
various attempts to correlate Kz with the Richardson 
number. Hunter and Spaulding applied Pritchard *s equation 
(62) for calculating vertical diffusion in the Providence 


River. 


2, ,, , „ „ ,-2 . (H-z) WH 

K = -r — (H - z ) (1 + B R. ) + a z ~ — — — 

z P ^ p H WT 


EXP (:|^) (1 t PpR.)-2 


(6,16) 


where 


, a , are adjustable constants R. = Richardson # 

P P P 1 

z = distance downward from surface H = total depth 

WH = wave height ^ horizontal velocity at 

WT = wave period 

WL = v7ave length 




The problem with Rg. 6.16 is that the physical density 
field and wind field must be measured, and this may indeed 
be difficult when trying to estiaate diffusion coefficients 
in offshore areas. In addition, the eguatioh has only been 
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SUMMARY OF VALUES 

OF VERTICAL DIFFUSION COEFFICIENT K- IN THE OCEAN 

C5) 

Note: Molecular diffusivity for heat; 1. 5 x lO’^ cm^/sec (at 20°C, 1 atm) 

Balt: 1.3 X 10”^ cm^/sec (at 20°C, 1 atm) 


Current or oceanic 
region 


Depth of 
layer (m) 


Vertical Diffusion 
Coefficient Ky 
(cm^/sec) 


Reference 


Philippine Trench 
Algerian Coast 
Mediterranean 


5000-9788 
0 - 20 
0- 28 


2. 0-3. 2 
35-40 
42 


Schmidt, 1917 
Schmidt, 1917 
Schmidt, 1917 


California Current 

0- 200 

30-40 

McEwen, 1919 

Caspian Sea 

0- 100 

1-3 

Stockman, 1936 

Barents Sea 


4-14 

Subov, 1938 

Bay of Biscay 

0- 100 

2-16 

Fjeldstad, 1933 

Equatorial Atlantic 
Ocean 

0- 50 

320 

Default, 1932 

Randes fjord 

0- 15 

0. 1-0. 4 

J acobsen, 1913 

Schultz Grund 

0- 25 

0. 04-0. 74 

Jacobsen, 1913 

Kuroshio 

0- 200 

30-80 

Sverdr up-Staff, 1942 

Kuroshio 

0- 400 

7-90 

Suda, 1936 

Southern Atlantic 
Ocean 

400- 1400 

5-10 

Defant, 1936 

Arctic Ocean 

200- 500 

20-5C 

Sverdrup, 1933 

Carribean Sea 

500- 700 

2. 8 

Seiwell, 1938 

South Atlantic Ocean 

3000-Bottom 

4 

Defant, 1936 

South Atlantic Ocean 

Near Bottom 

4 

Wattenberg. 1935 

West Atlantic Trough 
(50°S to 10°N) 

Near Bottom 

7-50 

Wust, 1955 


North Atlantic 
Indian Ocean 
Pacific Ocean 
Tidal Channe 1 

(Mersey estuary 
and Iri sh Seai 
Near Cape Kennedy, 

Florida 

Bikini Lagoon 

Coast of Denmark 
California Coast 


Near Bottom 

0 - 20 
(bottom) 

Surface Layer 
0- 50 

(bottom) 


2-40 

19 (in August) 
I. 3 (in Summer 
260 

0. 05 -T ■ " 

0 . 1-10 

15-180 

(at wind force 
3-4) 


Koezy, 1956 

Bowden, 1965 
(with R. from 
0. 1 to ^ 2. 0) 
Carter and Okubo, 
1965 

Munk, Ewing and 
Revelle. 1°4Q 
HarremoeSi 1 *“d7 
Foxworthy, Tib / 
and Barsom. lo66 
Stommel and 

W oodcock, 1951 


* As gi^^cn by Defant, 1961 
As given by Bowden, 1962 
As given by H ar r e moc s. 1967 
As given by Wicgel, 19o4 


TABLE 6.1 


Sunutiary of Values of Vertical Diffusion 
Coefficient in the Ocean . . , . . . 


- ij 
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SUMMARY OF FORMULAS 

ON 'CORRELATION OF VERTICAL DIFFUSION COEFFICIENT 

WITH RICHARDSON'S NUMBER R. (OR DENSITY GRADIENT e) 
(Reference S) ^ 

Note: K : K at =0, i, e., the neutral case 8 : proportionality constant varies 

^ ^ from case to case 


Rossby and Montgomery 
(1935)>t' 

= 

K 0 

Rossby and Montgomery 
(1935)* ■ 

''z: = 

K^od+SRi)-^ 

Holzman (1943)* 

(1 


Yamamoto (1959)* 

II 


Mamayev (1958)* 

II 

„ -8 a. 

Munk and Anderson 
(1948)** 

II 

N 

8 = 3.33 based upon data by Jacobsen (1913) 
and Taylor (1931) 

Harremoes (1988) 

> 

. ,--3 -2/3 2, 

5x10 X £ cm /sec 

note : £ in m ^ ; approximate experimental 
-<i -5-1 

range 5x10 <£ < 15x10 m 

Kolesnikov, et al 
(1961 )*** 

(1 

Q 2 

K, . + ^ in cm /sec 

r min c 

K , and 8 are empirically determined 

z min 



to be : 

K_ . = 12, 8 - 8, 3 X 10“^ (1958 and 

1 960 observations) 



K . = 2, e = 10.0 X 10'^ (1959 

2 min abse rvations ) 

Koh and Fan (1969) 


10 ^/G (K in cm^/sec; r in m"^^ 

^ - 7 -2-1 

4 X 10 ' < e ^ 1 0* m" 


■'•/is piyen by Oku bo ( 1 ) 

As pivfn by Buvvrlpn (I'^bii) 

The formulas prosonled in the translated version are aooarentlv erroneous , 


TABLE 6.2 


Summary of Formulas on Correlation of 
Vertical Diffusion Coefficient with the 
Richardson Number, (or density gradient) 
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verifi€>d for the James River. 

Csanady (33) remarks that the assumption that the 
vertical mean velocities and turbulence levels are small in 
comparison to those in the hori 2 ontal is true only when the 
wind is relatively light. He cites the existence of 
Langmuir circulations ( after Langmuir, (63) ) . When strong 
winds are present ( Langmuir circulations are invariably 
present for winds above 16 ft/sec), and in particular when 
the surface is strongly cooled, large scale horizontally 
aligned vortical flow structures, with axes parallel to the 
wind, sweep water downward from the surface to the first 
stable sheet (64,65). These are commonly known as 
windrows, being horizontally parallel rows of vortices 
having opposite rotational direction, with areas of 
upwelling and downwelling alternating between the rows. 
Fcam and debris accumulate at the lines of confluence of 
the Vortical rows. ( Figure 6.2 ) Their exact physical 
mechanism is not known (66), but whatever the mechanism, 
they rapidly distribute a ’’marked fluid” over the vertical 
range from the surface to the first stable sheet (33) . 

This mode of transport is not a turbulent mechanism, 
but is an advective flux, and cannot be described by a 
Richardson type vertical diffusivity coefficient 
formulation. If Langmuir circulations are present, then 
the effect is to rapidly distribute a marked fluid from the 
surface down to the first stable sheet. 

For more practical purposes of calculating sediment 
transport, a conservative estimate can be made, if no field 







data is available, that the water column is well mixed, and 
a constant vertical eddy diffusivity can be be estimated by 
Eq* 6.13d, or others given by Hunter and Spaulding (7) in 
Table 6.1 . If enough data is present, the formulation of 
Pritchard, Eg. 6,16, may be Used. 

When considering the turbulent transport of 
non-buoyant sediment particles, the specification of the 
vertical turbulence; level takes on new importance. 
Turbulence suspends the particles within the water column, 
just as dust particles are suspended by the " Brownian” 
random molecular motions. 

For sedimentary particles, the rate of turbulent 
diffusive transport of mass of concentration c in the 
fluid, may be written 


fK, 


^2 

0 C 




( 6 . 17 ) 


where beta is a constant on the order of unity, which is 
specific to a given constituent * s physical characteristics, 
such as diameter and roughness. Experiments by Vanoi (44), 
later confirmed by Brush (45) , and Milano (46) , indicate 
that beta has a value equal to one for sediment grain 

diameters less than 0.2 mm. for a dilute suspension less 

than 3 % by weight of sediment. 

In still, non-tur bulent water, if it were not for 
vertical turbulence, negatively buoyant particles would 
drop out of the water cloumn with a nearly constant rate 

equal to the asymptotic free fall velocity. Murray (67) 
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demonstrated that the vertical turbulence levels produce 
oscillatory vertical particle motion for entrained sediment 
particles which in effect retards the fall out rate from 
that predicted by simple asymptotic free fall theory. 

When trying to predict entrained sediment motion, 
WAPIC calculates a spatial-time averaged vertical diffusion 
velocity proportional to the product of the diffusion 
coefficient and the concentration gradient ( see Kg. 3.6 c 
) . 


w_ 


D 


K 

z 

°AVG 


'h c 
^z 


(6.18) 


A suspended particle is assumed to move with a vertical 
velocity which is the sum of the free fall settling 
velocity and the total vertical pseudo transport velocity. 
For purposes of illustration, assume that a channel 30 feet 
deep exisls with a uniform Xdirected mean velocity of 1.0 
ft /sec. 

The diffusion coefficient may be estimated using a 
iSanning coefficient of 0.020, and from Eg. 6.15-d, a K'z= 
0.09 ft2/sec. A particle of 0.000083 ft (0.0025 cm) mesh 
diameter with a specific densitY of 2.6 yields a settling 
velocity of 0.003 ft/sec from Rabey's modified settling law 
{ Eg. 4.13 ) Assuming a vertical mixing length of 3.0 feet 
( a very conservative estimate, which equals 30ft/10 HAPIC 
vertical grids) a diffusion velocity may be approximated as 


MAX 


_ 0.09 Ft /sec 
3 Ft . 


= 0.35 Ft /Sec 


(S.19) 
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Indeed, this value is ten times larger than the sediment 
free fall velocity. A large diffusion coefficient with a 
corresponding small vertical settling velocity may result 
in large excursions for a WAPIC particle leading to 
oscillations in the vertical concentration field if the 
simulation time step (or vertical accuracy number) is not 
sufficiently small. This effect will be discussed in a 
section of the following chapter. The choice of vertical 
diffusion coefficient will help to determine the fallout 


rate in the simulation of sediment transport . 
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VII. VEBIFICaTION RGAT-NSI ANALYTIC SOLUTIONS 
Before? WAPIC could be applied in the simulation of 
mass transport of sedimentary and pollutant materials in 
the coastal zone, some preliminary verification work was 
necessary. The purpose of the following studies was to 

test WAPIC and gain experience in the adaptation of WAPIC 
to various three dimensional mass transport situations. 
Five verification studies against analytic solutions to the 
mass transport equation were performed. The various 
operating parameters for each test case are summarized in 
Table 7.1 . 

A. Unbounded Unidirectional Turbulent Flow Field 
A spherically symmetric three dimensional Gaussian 
marker particle distribution was instantaneously released 
into a steady-uniform X-directed velocity field. Two 
thousand marker particles were placed, with a normal 
distribution in three space, within a 15x15x10 K^Y-Z 
Eulerian grid. A standard deviation to grid spacing ratio 
of 2.0 in the horizontal { X-Y ) was utilized in order to 
approximate, to within 6% , the correct diffusion velocity, 
thus insuring the propagation of the initially spherical 
marker particle distribution in time. Uniform and time 
independent eddy diffusion coefficients were utilized. 



BOUNDARY 

CONDITIONS 

GRID 

GRID 

EXP. 

TRANSL 

X Y Z 

X Y : 


NUMBER I INITIAL 

OF I STANDAW) 

PARTICLES I DEVIATION 

FT. 

a „ o 

xo ^yo zo 

X y z 


DIFFUSION I MISCELLANEOUS 

OR DISPERSION I OPERATING 

COEFT. I PARAMETERS 

FT^/SEC 


K K 

X y 


iO I 0.0 I Y y Y Y y Y I 15 


2000 I 1000 1000 1.5 50.0 50.0 10001 I Q * 10 


iO 0.0 Y Y Y Y Y Y 15 


2000 I 1000 1000 1.5 


n ■ 0.0001 sec 

50.0 50.0 0.0001 « 1.0 FT/SEC AT CENTER 

° OF CLOUD 


10.0 0.0 N N N N N N 20 


2816 240,0 240.0 3.0 


50.0 50.0 .0001 „ ,„8 CiQ 


1930 Secs. 


# PARTICLES/STEP = 88 


1.0 15 

5.0 0.0 Y N N Y N N 15 

0.0 15 


2000 I 60. 60. 2.0 


REFLECTION AT 
Y-Z BOUNDARIES 
W = 450 FT 
H * 30 FT 


0 1,0 

0 0.0 N N N y N N 

0 0.0 


3000 NA NA NA 


MANNING COEFF. “ 0.020 
0.0 7.7 H = 30 FT 

EQOA PARTICLE SETTLING RATES 

0.0005, 0.001, 0.005, 0.008 
0.01 FT/SEC 
REFLECTION AT SURFACE 
DEPOSITION AT BOTTOM 


NA = DOES NOT APPLY 
y «* YES N = NO 


n = SHEAR RATE (SEC ) 
= DEPTH W » WIDTH 
Q - MASS 


TABLE 7.1 Operating Parameters for Verification Against Analytic Soltitions 
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Comparison was made against a three dimensional 
analytic solution of the form, given by Quesada (68) - 


c(x,y,z,t) - 


+2Kt)V2(,2 

xo X yo y zo z 


X EXP (- 


(x - U t)' 
o 


Y 


2(0^ + 2K t) 2 (ct^ + 2K t) 2(a^ + 2K t) 

xo X yo y zo z 


•) (7.1) 


a = initial standard deviations 

xo,yo, *0 


u 


= uniform steady current velocity 
= distances from initial center of mass of cloud 


o 

x,y,z 
Q =5 mass 

K = asymptotic x,y,z turbulent diffusion coefficients 

x,y,z j. I u , 


Equation 7.1 states that the standard deviations of the 
cloud particle distribution in a homogeneous turbulent 
steady flow field grow in a " Fickian ” manner. 


a = (cf + 2K t) 
X xo X 


1/2 


(7.2) 


as the center of mass of the cloud moves downstream with 
the velocity n (33) . 

Figure 7.1 illustrates a two dimensional X-Y view oC 
the WAPIC particle distribution at time T=0, and 10000 
seconds. The axes of Figure 7.1 are in dimensionless X-Y 
grid units, where one grid unit ■= 500 length units. The 
center of mass of the cloud has translated a distance 
U^T = 10000 length units. Figure 7.2 details a direct 

quantitative comparison of WAPIC and the analytic solution 
for cloud concentration plotted as a function of downstream 
distance ( in grid units ) . The concentrations are 
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I plotted for the cloud, at various tines, along a line 

I • through the initial fflaximum concentration which is parallel 

I 

3 to the downstream x axis, WAPIC concentration predictions 

I are within an absolute value of 5 X of the analytic 

solution values. 

J : I 

I; An effort was made to illustrate the three dimensional 

j' ' 

i; diffusion process by producing two dimensional X-Y 

' ■ i i. ■ ' 

concentration contour plots at variouS| Z-levels, Figures 

“ 7,3 and 7.4 illustrate analytic and WAPIC contours of 

ir. ■ ■ . . ’ ■ r 

i concentration at times T=0,500d, and 10000 seconds, ( 



I Please ignore the difference in Z-levels in each: 

i| corresponding frame of Figs. 7,2 and 7. 3 ); This three 

dimensional unbounded diffusion study, verified WAPIC'S use 
of the translating and expanding grid options, and the 

ji j ' 

approximation to the one dimensional diffusion velocity. 


1 

I 


I 

li 
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J- 
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i 
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f 
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WBPIC PflRTlCLE DISTRIBUTIONS - INSTnNTRNEOUS GflUSS. RLSE. 


DIFFUSION IN UNIFORM X-FLOW VELOCITY=l. 0 
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FIGURE 7.x Tv;o Dimsnsional Particle Distributions at Time T - 0, lOQOO for Diffusion of an Initially 
Gaussian Release in a Unidirectional Steady Flow Field 





C/CMFIX ANALYTIC 

_p. 25 0. 50 0. 75 


WRPIC VERSUS fiNRLYTICfiL SOLUTION 




DIFFUSION - TNSTRNT. GRUSS. RLSE. IN UNIFORM X-FLOW 

CONCEN'TRBTIONS RLONG R LINE THROUGH THE INITIRL MRXIMUM 
CONCENTRRTION » PRRRLLEL TO THE DOWNSTRERM AXIS 
KX .KY ,KZ*=50. 0 .50, 0 lO. 0001 



X~D I STRNCE_ DOWNSTRERM ^ 

FIGURE 7.2 WAPIC Versus Analytic Solution, Diffusion an Initially Gaussian Release in a Unidirectional 
Flow Field, Concentrations Along a Line Through the Initial Maximum Value, Parallel to the 
Downstream Axis 
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CONCENTRRTION CONTOURS OKUBOS RNfiLYTICfiL SOLUTION 

DIFFUSION OF RN INSTRNTRNEG^U^^ RLSE. IN UNIFORM X-FLOW 

<TX CT2= 1000 .1000 >1. 5 DxDyDZ= 50. 0 .50. 0 .0, 001 
6X OY 02= 500 .500 .1. 5 


ZLEVEL=3.0 TSTEP= 5000.0 



U- 1.0 



ZLEVEL=4.0 TSTEP= 5000.0 



ZLEVEL-5.0 TSTEP= 5000.0 


ZLEVEL=6. 0 


5000. 0 


FIGURE 7.3 Co nc ent ration Contours, Okxibo's 2\nalytic Solution for Diffusion of an Instantaneous 
Gaussian Release in Uniform X-Flow 








CONCENTRHTION CONTOURS 


WRPIC 3D SOLUTION 


DIFFUSION CFRN INSTRNTRNEOUS GRUSS. RISE. IN UNIFORM X-FLOW 

“Tx lOCG .1000 »1. 5 DxDyDZ= 50* 0 .50. 0 .0. 001 

Ox 5y 5z=500 .500 . 1 . 5 

ZLEVEL=4.G TSTEP= 5000.0 ZLEVEL=5.0 TSTEP=‘ 5000.0 


U= 1. 0 



> cA 



ZLEVEL=6. 0 TSTEP= 5000.0 


ZLEVEL^r.O TSTEP= 5000.0 


FIGURE 7.4 Concentration Contours, WAPIC Solution to Diffusion of an Instantaneous Gaussian Release 
in Uniform X-Flow 
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B, Diffusion in a Dniformly Sheared Flow Field 
A three dimensional unidirectional velocity field with 
a uniform shear rate of du/dy = 0,0001 sec-1 , and a mean 
velocity of 1,0 ft/sec with a homogeneous turbulent field, 
was applied to study the shear effect on the diffusion and 
advection of an instantaneous Gaussian marker particle 
release. 

Comparison was made against the analytic solutions of 
Carter and okubo (69) and Okubo and Karweit (70) , with the 
modification to represent the diffusion of a finite 
Gaussian release in space as opposed to a point release. 
The analytic equation is of the form 


c(x,y,z,t) = 


0 

(2n)^'^(0^ + 2K t)^'^ (a^ + 2K + 2K t) 

xo X yo y zo z 


1/2 


(x - U t - ^ (fi y + Q zXt)^ 

EXP (- i ( ^ ^-X 5^ 

<°xo + ^3 ‘ > 


+ 2K t) 

yo y 



2 

z 



+ 2K t) 
z 



(7.3) 


where 


V, = 




K 


K 


2 ^ 
5 K 


)) 


(7.4) 
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The term ^ can be interpreted as the time at which the 
velocity shear begins to significantly affect the cloud's 
diffusion, or the time of the shear effect (70) . 

The shear effect can be seen in Figure 7.5, which 
illustrates a WAPIC two dimensional view of the particle 
distributions at time T-0 and 10000 seconds later. These 
plots show that the former Gaussian distribution is skewed 
to the side where the advection velocity is larger. The 
marker particle cloud disperses more ( or skews) with the 
levels of higher X-directed velocities. 

Figures 7.6 and 7.7, detail X-Y planar concentratioh 
contours at various Z-levels for the analytic and WAPIC 
solutions. Cloud concentrations are plotted at times T= 
0,5000 and 10000 seconds. The observation is made that the 
concentration contours are stretched in the direction of 
increasing velocity, thus enhancing diffusion in the 
lateral, Y, direction. The WAPIC predicted concentration 
value are within an absolute value of 5 % of the analytic 
values at time T= 10000. 


I 




CONCENTRRTION CONTOURS 


OKUBRS RNRLYTICRL SOLUTION 


DIFFUSICN - INSTRNT. GROSS. ROSE. IN UNIFORMLY SHERRED X-FLOW 


CTx ay 0*2= lOQO »1000 >1. 5 DxDyDZ= 50. 0 .50. 0 .0. 00 1 
5x 5y 500 .500 .1, 5 


ZLEVEL=F.O TSTEP^ 5000.0 



10.00 20 . 00 30 . 00 


Velocity « 1,0 0 < 


ZLEVEL=5.0 TSTEP= 5000.0 




10 . 00 20 . 00 30 . 00 


ZLEVEL=6.0 TSTEP= 5000.0 


ZLEVEL=7.0 TSTEP= 5000.0 



FIGURE 7,6 Concentration Contours, Okubo's Analytic Solution for Diffusion of an Initially Gaussxan 
Release in Uniformly Sheared X-Flow 
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C. Diffusion of Continuous Gaussian Releases 
A plume or continuous waste release in the long term 
transport phase may be simulated with WAPIC by releasing a 
Gaussian particle distribution at each time step (28,29). 
0kubo*s analytic solution for diffusion of a three 
dimensional instantaneous release in a uniform turbulent 
X-directed flow field, can be numerically integrated in 
time- to yield a solution for the diffusion of a succession 
of point or Gaussian releases. Since the center of the 
waste or particle release is considered to be fixed in 
space, the WAPIC Eulerian grid will also be fixed in all 
three dimensions in the non-expanding, non-translating 
mode. Particles reaching the boundaries will be removed 
from the computational field. 

An extremely large Eulerian grid of 45x20x15 X-Y-Z 
grid cells was required to perform a simulation of 1920 


seconds. 

with 

a mean 

velocity of 1, 

0 ft/sec and a X-Y grid 

size of 

120 

ft. 

Approximately 8 

8 marker particles were 

released 

per 

WAPIC 

time step. H 

orizontal standard 

deviation /grid 

spacing ratios of 

2.0 were utilized. 

Fig ures 

7.8 

and 7. 

9 illustrate th 

e concentration co n t o ur 


plots for the analytic and WAPIC plume simulation. 

WAPlc shows correct plume concentration contour trend 
behavior, but shows irregular contour lines with 
considerably less diffusion along the outer fringes, 
relative to the analytic solution. Concentrations near the 
center of the WAPIC release were 3-4 times greater than 
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that found in the analytic solution. The reason for this 
poor behavior is that when using only 88 particles to 
represent each release, the particles tend to be cluroped 
near the release point and not ” Normally ‘' distributed in 
space. This causes high concentration gradients and poor 
spread along the fringes of the release cloud where there 
are very few particles. As explained in Chapter V., this 
problem could have been overcome through the use of 
approximately 1000 or more marker particles per time step. 
But due to run and computer storage limitations, this would 
require at least 512k bytes of computer storage, to achieve 
an accurate comparison with the analytic solution. The 
inference can be made that the WAPIC method is not a 
practical method for simulating a continuous waste release 
unless the user has an extremely larga, fast computer. 
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FIGURE 7.8 Concentration Contours, Analytic Solution for Diffusion of Continuous Gaussian Release in 
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D. Diffusion in Unifora Bounded Channel Flow 
This simulation was performed to study W^PIC'S ability 
to model diffusion which is horizontally bounded^ ie. 
reflecting boundaries. Analytic solutions are detailed in 
the works Of . Gsanady (33) , Wnek and Fochtaan (71), and 
Cleary and Adrian (72) . Csanady and wnek utilize the 
method of images to approrinate diffusion behavior at 
reflecting boundaries, while Cleary and Adrian use the 
integral transform approach, whose solution may be written 
as 


C(x,y,z,t) 


X - u t) 

gSXp 3^ 


WH (4TTK^t) 


1/2 


X (1 + 2 S exp (- K t) cos (P„z) ) 

N z N 

n=l 


(7.5) 


X (1 + 2 S exp (- K t) cos (Bj^y)) 
n=l ^ 


where - 


isfrr MTT 


H = depth, W = width 


A Channel 30 ft deep and ^i50 ft, wide with a uniform 
steady current field of 1.0 ft/sec was chosen as the 
experimental "apparatus" , The WAPIC Rulerian grid system 
was allowed to translate and expand only in the X-direction 
while the outer edges of the 1 and Z grid dimensions became 
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the channel boundaries. 

Cleary and Adrian's solution technique applies to 
diffusion Erom a point release. In order to simulate a 
Gaussian release at time T=0. the analytic solution had to 
be advanced in time 72 seconds, assuming a uniform 
turbulent flow field in the horizontal and vertical. This 
estimate follows from the Ficklan turbulent diffusion 
relation 

= (2k^ (7.6) 

for an instantaneous point release to grow to a given 
standard deviation. 

Concentration contours for VAPXC and the analytic 
solution at different vertical levels for times T=0 and 
T=720 seconds later are given in Figures 7.10 and 7.11, 
respectively. The observation is made that the effect of 
the lateral boundary is to equivalence values of 
concentration in the lateral, f, direction as the release 
moves downstream. For a line drawn through the center of 
the channel in the X-direction, concentrations agree within 
5% . Along the walls the deviations are greater, in the 
5-10% range. This large deviation can be attributed to the 
initial placement of the mar)cer particles with reference to 
the channel center. The analytic solution has the point 
release placed exactly in the channel center, with a 
radially uniform concentration distribution^ which results 
in symmetric reflective behavior at time T-720 seconds. 
Because WAPIC’S initial particle distribution is not 
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perfectly symmetric around th® channel center, perfect 
lateral symmetry of concentration reflection is not 
possible . 
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FIGURE 7.10 Concentration Contours, WAPic Solution for Diffusion of an 
Instantaneous Gaussian Release in Uniform Bounded Channel 
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E. Simulation of Sediment Transport in Stream Flow 

Since WAPIC will be utilized to simulate the long term 
transport and fallout of suspended sediment in the coastal 
and estuarine environment, some preliminary testing of the 
pseudo total velocity method for treating turbulent 
settling was necessary. The decision was made to simulate 
the case of sediment fallout from an initially vertically 
uniform one dimensional concentration distribution in 
turbulent steady unidirectional channel flow. In 1943, 
Dobbins (73) derived an analytic formulation which 
successfully explained the observed vertical sediment 
distribution in steady one dimensional stream flow, in 
which the dilute suspended sediment was in equilibrium with 
the bed material. He also developed an analytic equation 
expressing concentration as a function of the vertical and 
horizontal downstream distances for an initially vertically 
uniform sediment distribution falling to zero concentration 
with the no scour bottom boundary condition. Sayre and 
Chang (53) integrated this equation over the vertical to 
obtain the vertically averaged concentration, or mass, 
remaining in the water column as a function of downstream 
distance. 

A review of the equations and simplifying assumptions 
used by Dobbins to produce this solution would be 
instructive. Consider a turbulent unidirectional ( X-Z ) 
sediment laden flow field of singular settling velocity 
bounded by a uniform free surface and hot tom boundary 



W = PARTICLE SETTLING VELOCITY 
S 


Dobbins assumes that the X-velocity distribution in the z direction, though 
logarithmic and nearly parabolic, is approximately uniform across the verti- 
cal, and the vertical diffusion coefficient is uniform with depth and can 
be related to the bottom shear stress or shear velocity by 

K = kHU* (7.7) 

Z : 

U* = g^ 

C = 1.49H^'^®/C (7.9) 

h m 

where : ^ 

C = Manning Coefficient 
m 

C^ = Chezy Coefficient 

k = Von Karmann Constant = 0,40 
H = Hydraulic Radius or Depth 
U* = Shear Velocity 
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Equation 7.7 relates the matheraatical expression for 
vertical eddy diffusivity for a parabolic velocity 
distribution . Equations 7.8 and 7.9 are the Chezy 

expressions for relating shear velocity to a bottom 
friction coefficient, the Manning friction factor. 

The two dimensional mass transport formulation may be 


written as 




z^2 


.2 ^2 
1L2. + K hs. 




( 7 . 10 ) 


Assuming that the longitudinal diffusion term is small in 
comparison with the vertical diffusion 


- + w ^ ° 


( 7 . 11 ) 


Since the fluid is moving with a uniform velocity, U , the 
Lagrangian viewpoint of an observer moving with the 
fluid-sediment parcel at velocity 0 is taken. The 
assumption is made that 


Since 


X = ut 


^x = u^t 


and Eg. 7.11 may be rewritten as 


^ c c ^ c 


( 7 . 12 ) 


Equation 7 , 12 describes the change of concentration of a 
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turbulent fluid parcel containing a discrete settling 
velocity, Hs, as a function of time* This linear partial 
differential equation was solved for the boundary 
conditions of no net sediment transport across the free 
surface and for no scour at the bottom for any given 
initial vertical sediment distribution. The solution for 
particle concentration as a function of depth and time ( or 
distance downstream) for the case of dropout to zero 
concentration from an initially vertically uniform 
distribution is 


/ 4.V ^ -0.5W 2(K ) 

c (z,t) = C e s z 

o 


-1 


„ 2 K + 0.25 W ^ K ”^)t , 

® 2o e n z s z 4 

s. EL ” 

n=l 2 ^ ^ ^ 

(O n+W^/4K^)(«.n t^) 

4k z 

z 


(7.13) 


where are the real successive positive roots of the 

transcendental equation 


2 cot (Hey) = Ho/CW H/2K ) - (W H/2K )/(HCk') 

s z s z 


(7.14) 


and ij/ is defined as 


w 

g 

t = cos{cy z) + (-rr ^) sin (a z) 

n n 2K cy n 

z n 


(7.15) 


Sayre and Chang integrated Eg. 7.10 over the vertical 
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to yiolfl the ratio of the average concentration in the 
vertical to the initial average concentration. 


K xU* 

(-1)""^ a ^ (-9p^ + 

n exp n (6HU_) 


(c/c ) = 72p^e^^ £ 

(9p^ + a_^ + 6p) (9p^ + 


Where Ut, the shear velocity is approx in a ted by Eg. 7.8, 
and ^ f the vertical dimensionless fall velocity is 

defined as 


MG# m POiKp. 


Equations 7.13, 7.14, 7.15 and 7.16 were computer coded 
to solve for the concentration distributions as a function 
of the vertical coordinate and tine, and for the fallout 
rate respectively. 

Humerical experiments were performed by WAPIC in order 
to simulate solutions to Egs. 7.13 and 7.16 . A channel 30 
feet deep was divided into 15 vertical grid cells with a 
3x3 X-y Eulerian grid matrix. Three thousand marker 
particles were horizontally centered in the exact middle of 
the X-y matrix, and were equally spaced from the free 
surface to the bottom in order to develop an initially 
uniform concentration disttibution for one coluran of 
vertical cells. A unidirectional X-flow field of 1,0 
ft/sec was specified throughout the Eulerian grid. A 




constant, uniform vertical diffusion coefficient was 
calculated by the relationships in Egs. 7.7, 7,8, and 7.9, 
with a Hanning friction factor of 0.020 yielding a 
turbulent diffusivity, Kz, of 0.09 ft^/sec. The X-Y 
turbulent diffusion coefficients, Kx and Ky, were set 
identically to zero. The Bulerian grid was allowed to 
translate only in the X-direction with the mean velocity. 
The free surface was given a reflection boundary condition: 



and the marker particles were removed from the | 

s 

computational field as they advected across the bottom. | 

For the initial simulation a settling velocity of | 

a 

I 

0.005 ft/sec was assigned to each particle. Figure 7.13 | 

i) 

demonstrates the effect of MAPIC time step on the | 

. I 

prediction of mass fallout rate as a function of distance f- 

downstream, for T = 60, 30 and 20 seconds. The observation | 

is made, that as time step decreases, the WAPIC solution | 

technigue approaches the analytic solution. At time steps p 

of 60 and 30 seconds the concentration field was found to | 

■ - I 

oscillate about the analytic solution. f 

These results, plus the fact that WAPIC employs a | 

■ ■ 

space centered finite difference approximation to define P 

, 'I 

the total pseudo velocity at the cell faces, suggested that | 

a vertical numerical stability requirement of the following ii 

' ■ v. 

form should be utilized for marker particle movement (33) t; 
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COMPARISON OF NARIC AND ANALYTIC MODEL FOR PREDICTION OF 
CONCENTRATION AVERAGED OVER THE VERTICAL VS DISTANCE DONNSTRERM 
FOR WHPIC TIME STEPS OF 60 j30 AND 20 SECONDS 
SETTLINC VELOCITY^O. 005 FT/SEC, VERT. DIFF. COEFF. =0. 09 FT2/SEC 
PARTICLE DEPOSITION IN UNIFORM TURBULENT X-FLOW , U=l. 0 FT/SEC 



FIGURE 7.13 Comparison of WAPIC and Analytic Solution, Effect of WAPIC Time Step on the Prediction of 
Fallout Rate as a Function of Dovvnstream Distance for a Singluar Particle Settling Rate 
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Vertical ^ 
Accuracy 


W At 



^ z 


K A t ^ 

+ £0.5 




( 7 . 19 ) 


Equation 7,19, which is a '• Courant" dynamic stability 
requirement or time step restriction, states that for 
stability of the solution technique, a marker particle 
cannot be moved a distance greater than one-half a grid 
length per time step. Hs, can be visualized as a constant 
vertical advective velocity. 

Table 7.2 details values of the vertical accuracy 
number ( Equation 7.19) versus the WAPIC time step for 
values of Kz=0.09 ft2/sec, Ws=0.005 ft/sec, and vertical 
grid size, AZ =2.0 ft. 

Time Step -sec. Accuracy number 

60 
30 
20 

TABLE 7.2 


1.50 
0.7 5 
0.50 


Figure 7.14 details the development of an instability 
in the MAPTC predicted concentration field as a function of 
time or distance downstream. The instability starts at the 
free surface where the diffusion gradients first initiate, 
and then propagates through the vertical concentration 
field, until the whole concentration field oscillates at 
further distances downstream. 
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riguce 7.15 shows the predicted vertical concentration 
field, with a WAPIC time step equal to 20 seconds, or 
accuracy number of 0.5. The concentration oscillation has 
been effectively damped out and the technique predicts the 
concentration field within an absolute accuracy of 4.0 % . 

Utilizing a time seep of 20 seconds, WAPIC simulated 
the fallout rates of five particle settling velocities of 
0.0005,0.001,0.005,0.008, and 0.01 ft/sec, with the same 
operating parameters used as above. Figure 7.16 details 
the fallout rates as a function of distance downstream. At 
the lower settling rates of 0.0005 and 0.001 ft/sec, the 
plots demonstrate that the fallout rate is essentially the 
settling velocity of the particles within a one hour time 
span from time T=0. This behavior indicates, for this 
case, that the smaller free fall settling velocities are 
rate limiting, not allowing the dominant diffusion process 
to come into action. As the free fall particle settling 
velocity increases relative to the particle diffusion 
velocities, for the diffusion value shown, the deposition 
rate decreases because concentration gradients develop near 
the free surface. The diffusion velocity acts to move the 
marker particles back into regions of depleted 
concentration which decreases the fallout rate. 
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Damping of Ihe Instability in the WAPIC Concentration Field by Use of the Accuracy Number 
Requirement 
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VIII. SIMULATION OF LONG TERM SEDIMENT TRANSPORT 
A. Objectives and Background 

As explained in Chapters IX, and III, WAPIC was 
designed to simulate the three dimensional long term 
turbulent diffusion and advection of dilute, neutrally and 
non-neutrally buoyant sedimentary and pollutant 
particulates in coastal and estuarine waters. The method 
was utilized to simulate the long term transport of a 
dilute suspended sediment cloud resulting from a barged 
waste release at Brown’s Ledge in Rhode Island Sound. This 
area is currently being considered by the U. S. Array Corps 
of Engineers as a dump site for Fall River dredged spoil 
disposal (1). The dump site is a square nautical mile 
area, the center of which is located 41® 18.3*N, 71® 04.1^W 
in an area known as Brown's Ledge. ( Figure 8.1) From 
Figure 8.2, which details the bathymetry at Brown's Ledge, 
it can be seen that the area surrounding the dump site, is 
relatively shallow water ranging from 60 - 120 feet , the 
dump site being 100 to 110 feet deep. 

Since this area is relatively shallo", the hypothesis 
is made that the short term transport phase of the dredged 
waste release would be entirely dominated by the 
entrainment phase (6). In other words, practically all of 
the sediment would descend directly to the bottom, with 
very little loss of material to the ambient current system, 
Gordon (74) in an observational study of dredged waste 
disposal in Puget Sound, found that the material left the 
barge with a "well definad mass falling with a maximum 


FIGURE 8,1 Geographic Location of Browns Ledge and the proposed One Square 
Nautical Mile Dump Site 



FIGURE 8.2 Depth Contours at Browns Ledge 
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velocity Of approximately 6 ft/sec. <« Using this value as a 
representative estimate of the descent velocity of a waste 
cloud, at Brown's Ledge, the cloud would be expected to 
reach the bottom well within 60 seconds. Any sediment 
entrained in the water column would be due to the ambient 
turbulence, and the surge of the cloud as it lands on the 
bottom . 

A very conservative engineering estimate can be made 
that no more than five percent of the instantaneous dump of 
bottom sediment would remain within the water column to be 
advected and diffused by the tidal and wind driven current 
fields. Since the average sea scow carries approximately 
3000 cubic yards (2) , and the typical specific density of 
Pall River sediment being 2.6, then approximately 10® grams 
(110 tons) of sediment would be entrained within the water 
column. If this cloud initially possessed a cylindrical 
shape 200 ft. in diameter and a vertical length of 35 feet, 
and the sediment were distributed evenly throughout at time 
T = 0, then a mean concentration of 3000 mg/liter (0.187 

lb/ft3) would be apparent. This value is within the limits 
of three percent by weight for dilute settling and 
transport of suspended sediment noted by Vanoi (44) . 

This suspended sediment will advect and diffuse with 
the tidal and wind driven currents in Rhode Island Sound. 
As described by Saila et . al . (2) , two significant field 

studies of the hydrodynamics of Rhode Island Sound have 
been accomplished. The non-tidal wind driven circulation 
was investigated by Cook. (75) in 1966, using wind drift 
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bottles and sea bed drifters to determine upper and lower 
layer net drift. His generalized net drift data for Rhode 
Island Sound are detailed in Table 8.1. 



SURFACE 

BOTTOM 

SEASON 

DRIFT 

DRIFT 

Spring 

N and E 

NW 

Summer 

N 

NW 

Autumn 

S 

N 

Winter 

. S 

N 

SPEED 

1.1 to 7.6 

0.05 to ; 

RANGE 

n.m./day 

n.m./day 


Table 8.1 Cook's 1967 Study of Non-Tidal Net Drift in Rhode Island 
Sound 


Shouting (76) in 1967 utilized an array of current 
meters to perform a 13 day study of the horizontal current 
patterns within a square mile area in Rhode Island Sound. 
Shonting found that tidally driven inertial flows combined 
with tidal motions dominate the current regimes and that 
mean speeds in the upper layer of current motion did not 
exceed Q. 50 knots, and that turbulent eddies were present 
having scales less than a kilometer. 

Hurlburt and Spaulding (8) applied a two dimensional 
tidally driven hydrodynamic model to predict current 
patterns and magnitudes in Rhode Island Sound - Block 
Island sound - Buzzards Bay Complex. An interesting result 
of this modeling effort was the indication of extremely 
weak tidal currents in Rhode Island Sound, especially in 
the region of the proposed dump site. Figures 8.3 A, B, C, 
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TIDAL CURRENTS IN KNOTS ZERO HOURS 
AFTER HIGH WATER AT NEWPORT, RI 


FIGURE 8. 3 A Two Dimensional Vertically Averaged Velocity Vectors in Rhode 
Island Sound, 0.0 Hours After High Water at Newport, R.I, 
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FIGURE 8.3B Two Dimensional Vertically Averaged velocity Vectors in Rhode 
Island Sound, 4.0 Hours After High Water at Newport, R.I. 
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AFTER HIGH WATER AT NEWPORT, RI 


FIGURE 8.3C Two Dimensional Vertically Averaged Velocity Vectors in Rhode 
Island Sound, 8.0 Hours After High Water at Newport, R.I. 



TIDAL CURRENTS IN KNOTS TWELVE HOURS 
AFTER HIGH WATER AT NEWPORT, RI 


FIGURE 8. 3D . Two Dimensional Vertically Averaged Velocity Vectors, in Rhode 
Island Sound, 12,0 Hours After High Water at Newport, R.I. 
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D detail tidal currents in I^hode Island Sound, referenced 
to slack water at Newport, Rhode Island. These figures 
show that the vertically averaged tidal currents range from 
0.1 to 0.4 knots during the twelve hour tidal cycle. 

Utilizing data from this modjl, an average net drift 
in the proposed disposal area was found to be approximately 
0.023 knots or 0.55 n.ra./day, in a direction measured 110 
CW. from true north {almost ESE) . The magnitude of this 
value is within range of values found by Saila et . al . (2) 
in July, 1970 of 0.03 - 0.08 knots net drift with a general 
southeast to southerly direction. 

B. Procedure of Application 

The transformed version of WAPIC was coupled to 
Spaulding and Hurlburt's (8) tidally driven two dimensional 
vertically averaged hydrodynamic model through the use of a 
two dimensional velocity processor. As explained 
previously, the velocity processor was required to 
interpolate advective velocities from the square nautical 
mile meshes to the initially smaller, translating and 
expanding three dimensional WAPIC Eulerian velocity grid. 

In order to adjust the hydrodynamic velocities with 
experimentally reported values for net drift at Brown's 
Ledge, a correction procedure was employed (77). The 
procedure was to subtract from the hydrodynamic velocity 
vectors, in the immediate Browujs Ledge area, the net tidal 
drift (calculated from the 12 hour average) at the proposed 
dump site. Subsequently, the net drift velocity vector was 
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added -to the new velocity components with a correction for 
the variation of bottom depth relative to the depth at the 
area of the experimentally determined net (wind and tidal) 
drift. The overall effect was to subtract out the model 
net tidal drift, and add the total net wind and tidal drift 
(based on field data) in the area to produce a velocity 
field yielding a net drift in the direction desired. The 
total net drift chosen for this particular set of 
sinnlations was 0,7 n.m. day in the ENE direction (77). 

Seven particular simulations were performed to 
illustrate the diffusive and advective patterns of mass 
transport in the area. The first three cases utilized 
marker particles with neutral buoyancy and constant 
horizontal diffusion coefficients of 50.0 and 100.0 ft^/sec 
for the first two and scale dependent horizontal diffusion 
coefficients for the third ( Table 8.2). Four additional 
simulations of vertical fallout of sediment were made by 
assigning settling velocities to the individual Lagrangian 
marker particles. Various mean sediment diameters and a 
specific density of 2.6 for Fall River sediment (1) were 
assigned to each marker to produce free fall settling 
velocities with the use of Eg. 4.13. 

All four of the marker particle settling velocity 
studies utilized constant horizontal diffusion coefficients 
of 100,0 ft2/sec. These four simulations were divided into 
two sets of two simulations. The differing variable in 
each simulation of both sets was the magnitude of the 
vertical diffusion which was 0.0001 and 0.01 ft2/sec 



X-Y GRID 
LENGTHS 
FT. 

DIFFUSION 

COEFFTCIENTS 

ft'^/sec 

INITIAL X,Y,T1 
STANDARD 
DEVIATIONS 
FT. 

# 

OF 

MARKER 

PART. 

PARTICLE 

DIAMETERS 

FT. 


X, Y 

K , K, 
X Y 


a 

X 

O 

y 



N.B. = NEUTRALLY 
BUOYANT 

1. 

400 

50. 

.0001 

800. 

800. 

ATI 

3000 

N.B. 

2. 

200. 

100. 

.0001 

400. 

400. 

ATI 

3000 

N.B. 

3. 

200. 

1 

Scale 

Dependent 

.0001 

400. 

400. 

ATI 

3000 

N.B. 

4. 

200. 

100. 

.0001 

400. 

400. 

AT] 

2 500. 

0.000083 FT 
(0.0025 cm) 

5. 

200. 

100. 

.01 

400. 

400. 

ATI 

2 500. 

0.00083 FT 

6 . 

200. 

100. 

.0001 

400. 

400. 

ATI 

2500. 

0.00005, 0.000083 
0.00016, 0.0002 
0.00028 FT 

7. 

200. 

100. 

.01 

400. 

400. 

ATI 

2 500. 

0.00005, 0.000083 
0.00016, 0.0002 
0.00028 FT 


TABLE 8.2 Operating Parameters for Browns Ledge Simulation of Long Term Transport 
of Sediment 
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respectively. The first set of simulations considered the 
transport of a single particle diameter of 0.000083 ft 
(0.002S cm), which yields a settling velocity of 0.003 
ft/soc by Eg. Q.13. The second set simulates the transport 
of five different particle diameters of 0.00005, 0,000083, 
0.00016, 0.00020, and 0.00028 ft. (0.0015, 0.0025, 0,0040, 

0,0060, and 0,0085 cm.) which yielded settling velocities 
of 0.001, 0.003, 0.008, 0.117, and 0.026 ft/sec 
respectively. 

The simulations of a single and five different 
particle fall velocities demonstrated the effects of 
choosing a single average particle diameter or fall 
velocity to represent the turbid cloud as opposed to a 
spectrum of particle diameters. The variation of the 
vertical diffusion coefficient was used to determine the 
effect of the vertical turbulence scale versus the particle 
settling rate on the fall velocity of the total particle 
cloud. 

For each of the seven simulations, the initial 
particle distribution was three dimensionally Gaussian in a 
15x15x10 X-Y-lfl Eulerian grid. The particle cloud was 
placed directly above the shelf at Brown’s Ledge, gust 
slightly to the northeast of the proposed disposal site. 
This location was chosen for simulation, because of the 
interesting variation in bottom topography ( Figure 8.2) 
and the economy of simulating particle fallout in 60-70 
foot depths as Opposed to 100 - 120 foot depths at the 
disposal site . 
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The center of mass of the particle cloud vas placed 

approximately 20 feet deep, the cloud possessing a vertical 

range of fifteen feet. All four settling velocity cases 
also utilized a 0,7 n.m./day net wind drift in the ENE 
direction. The negatively buoyant particle simulations 

were run until more than 95 per cent of the marker 

particles had fallen to the bottom. 


C. Results and Graphics of Simulations 

1. Kx, Ky= 50 ft^/sec. Neutrally Buoyant. 
This neutrally buoyant simulation was run for 


approximately 

15 

hours 

during 

wh ich 

an initial maximum 

concentration 

of 

3000 

mg/liter 

decayed to 1-5 mg/liter 

levels for 

constant horizontal 

and 

vertical diffusioti 

coe f f icients 

of 50 

. 0, and 

0.0001 

f tz/sec 

respectively . 


Figure 8.4 A details the translation and expansion of 
an initially square nautical mile WAPIC grid in time 
through the Sound waters in approximate three hour time 
steps. Figures 8.4 B, C, D, E, F, G illustrate two 

dimensional X-Y concentration contours at vertical levels 
of 4, 5, 6 and 7 at different times in the 15 hour 
simulation. Concentrations were plotted at these levels 
because they contaiwed the maximum particle densities. By 
time T = 12420 sec.. Fig. 8.4 C, the original HA PIC grid 
cell length had expanded from 400 to 781.2 feet. The 
skewness in the cloud was evident as the cloud grew large 
enough to feel the shear in the advective tidal velocities. 
By time T = 41,100 secs.. Fig. 8.4 G, the WAPIC grid had 
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translated and expanded to the lower edge of the adwective 
velocity field and had stopped expanding and translating, 
while the particle cloud was no longer located in the 
center of the Euleriaii WAPIC grid , 

2. Kx , 1^, = 100.0 f t^/sec. Neutrally Buoyant 
This second neutrally buoyant marker particle 

simulation utilized horizontal and vertical diffusion 
coefficients of 100.0 and 0.0001 ftz/sec respectively, and 
an initial X-Y grid size of 200 feet. Since the horizontal 
diffusion coefficients had doubled, the background levels 
were reached in less than half the time of the previous 
simulation. Figure 8.5 A details the two dimensional 
translation and expansion of the horizontal WAPIC grid at 
times T = 0, 1, 3, 5, 6, and 7 hours after slack tide. 

This plot demonstrated a net drift of the marker particle 
cloud in an ENE direction. With the original area expanding 
to more than 36 times its original size. Figures 8.5 B, C, 
D detail two dimensional concentration contours at various 
simulation times. After six hours, ( concentration contour 
plots are not given at T = 21600 sec. ) , the maximum 

concentration level had dropped to the 1-5 rag/liter level, 
as compared to fifteen hours for the previous simulation. 

3. Scale Dependent, Neutrally Buoyant 

A simulation was done to test WAPIG'S horizontal scale 
dependent diffusion coefficient calculation procedure. The 
horizontal diffusion coefficients were calculated by 
Leendertse's formulation. Eg. 6.13 a ^ This equation 
approximates the total diffusion coefficient as the sum of 
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Elder's diffusion Coefficient plus scale dependent 
diffusion due to the wind. At each tine step, WAPIC 

calculated the horizontal standard deviations and an 
updated by Eqs. 6,13 e, f. 

An initial grid size of 200 feet was utilized with an 
initial maximum concentration of 1500 mg/liter. After 
approximately 5.0 hours of simulation, concentration levels 
had fallen to background values of 1-5 mg/liter. Figure 
8.6 A details the horizontal grid expansions and 
translations at T = 0, 1.0, 1.8, 2.6, 3.6, and 4.6 hours. 
Again, the net drift of the particle cloud was towards the 
ENE, as would be expected. 

Figures 8,6 B and G show two dimensional concentration 
contours for an X-Y plane through the maximum concentration 
level at various times. The maximum concentration stays at 
the same vertical level. The values of the X-Y grid 
lengths, latitude and longitude of the WAPIC grid origin 
for each time level of concentration plot are given at the 
right hand side of Figs. 8.6 B, C in descending time order. 
In ccraparison with the two previous simulations, the 
utilization of the scale dependent horizontal diffusion 
coefficient approximation, accelerated the transport of the 
cloud concentration to background levels of entrained 
sediment. 

4. Single settling Rate, Kg =0.0001 ft^/sec 

The fourth simulation utilized 2500 markers, each with 
a free fall velocity of 0.003 ft/sec in a tidal velocity 
field, with Fickian horizontal and vertical diffusion 
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coefficients of 100.0 and 0.0001 ft^/sec respectively, 
figure 8.7 A details the expansion and translation of the 
horizontal grid at tines T = 0, 3, and '5 hours after 

high water. Figures 8.7 B, C, D illustrate two dimensional 
concentration contours for an Jr'i plane (for a given 
vertical level) passing through the maximum concentration 
at various times. These figures illustrate the descent of 
the regions of maximum cloud concentration with increasing 
vertical levels. Figures 8.7 E, F, G, H demonstrate the 

physical vertical descent of the marker particle cloud to 
the bottom. These are X-Z, and Y-Z planar views of the 
marker particle cloud in space. The axes are in units of 
feet ( depth ) and dimensionless horizontal grid distance. 
The intersection of a plane passing through the cloud's 
center of mass and intersecting the bottom, in either the 
X-Y direction, yields a bottom topographic line for 
referencing the marker particles. One out of ten of the 
marker particles were plotted. 

These figures show the physical descent and horizontal 
expansion of the marker particle cloud through time. 
Figure 8.7 H displays the final descent of the cloud at T = 
5.19 hours to the bottom. Particles below the bottom 
topographic line were not below ground, but were just below 
the vertical level of the intersection of the vertical 
plane and the bottom. 

Figure 8.7 I details the vertical depth and horizontal 
translation of the center of mass of the particle cloud 
from T = 0, and the percent mass deposited of the marker 


144 


particle cloud as a function of time. The deposition rate 
was nearly constant from time T = 3.0 to 5,0 hours, then 
decreased when the particles began to reach the bottom and 
were diffused and advected at a nearly constant vertical 
level. After 6,5 hours, more than 95 pecent had been 
deposited. The X-Y coordinates of the center of mass of 
the cloud translated in an ENE direction. The slopes of 
the first two plots of Fig. 8,7 I record the vertical and 
horizontal velocities of the center of mass while the third 
slope was the actual deposition rate, 

5. Single Settling Rate, Kg = 0.01 ft^/sec 
The fifth simulation was exactly the same as the 
fourth except a vertical diffusion coefficient 100 times 
larger, K 25 = 0.01 ft^/sec, was utilized. In Fig. 8.8 A, 

the grid translated towards the East, with a net 
displacement of the grid origin in the ESE direction after 
6,5 hours of simulation. Figure 8.8 B details horizontal 
plots of the Concentration contours, and 8.8 C, D, E, F 
show the vertical descent of the marXer particle cloud. 
Hithin a time span of four to six hours following the 

i I 

release. Figure 8.8 G, a lower particle deposition rate was 
evident as compared to that of Fig. 8.7 I. The plots of 
the vertical positions of the cOnter of mass for the above 
two figures were almost identical. The initial effect of 
the higher diffusion coefficient in case 5 was to increase 
the deposition rate (more vertical spread), but as the 
center of mass of the cloud reached the bottom, the effect 
was to retard the fallout rate from that found in case 4, 
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I nterestinqly , both caso 4 and 5 reached the 95 percent 
mass deposited level at the same time. 

6. Five Settling Rates, Kz - 0.0001 ft*/sec 

The sixth simulation utilized five different marker 
particle settling rates (i.e., diameters) with constant 
horizontal and vertical diffusion coefficients of 100.0 and 
0.0001 ft^/sec respectively. The 2500 marker particles 
were evenly divided into groups of 500 for each particle 
diameter. All particles were assumed to possess the same 
specific density of 2.6. 

Figure 8.9 A shows the horizontal translation and 
expansion of the WAPIC grid at times T = 0, 2, 4 hours. 
Figures 8.9 B, C, show the concentration contours for the 
maximum cloud concentrations at various times# T =0, 1.12, 
1.81, 2.16 hours. Figures 8.9 D, E, F again illustrate the 
vertical descent of the marker particle cloud, to the 
bottom. Figure 8.9 G in comparison with Fig. 8.7 F, both 

at the approximate tine of T ,= 2.2 hours, demonstrated a 
much greater vertical spread in particle position due to 
the fact that there were five particle settling rates, with 
three velocities greater than the settling rate in case 4. 
In case 6, f"'ie vertical particle movements were dominated 
by the magnitudes of the settling velocities which 
overshadowed the diffusive vertical flux. 

In Fig. 8.9 H, the particle cloud deposited its mass 
in less than 2.5 hours and moved no further than 0.7 
nautical miles to the enE of the original position. Figure 
8.9 H further demonstrated that the settling velocity of 
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the center of mass of the cloud varied sliqhtly with time. 

7 i Five Settling Rates, Kg - 0,01 ft^/sec 
The seventh, and final simulation of particulate 
transport at Brown's ledge was the same as the sixth except 
that a constant vertical diffusion coefficient of 0,01 
ft^/sec was , utilized. The purpose of using this larger 
coefficient was to determine the effect of the magnitude of 
the vertical turbulence scale relative to the settling 
cates on the overall fallout rate of the cloud. Figures 
8.10 A, B, C, D, E, F, G relative to Figs. 8.9 of case 6, 
demonstrate that there was very little difference between 
the two simulations as far as concentration, center of 
mass, and fallout rates were concerned. This behavior was 
due to the fact that the magnitudes of sixty percent of the 
particle settling velocities were greater than the 
diffusion velocities of the particles. In other words, for 
cases six and seven, the free fall velocities dominated the 
vertical scale of notion. These simulations demonstrate 
that adequate resolution of particulate fallout from a 
suspended sedinent qloud requires not only knowledge of the 
vertical scale of lurbulence but also the distribution of 
sediment dianeters and densities (fall velocities) . 


COUPLING WRPIC TO 2-D VERTICRLLY RVG. HYDRODYNBMIC TIDRL MODEL 
TRRNSLRTION RND EXPRNS ION OF 15X15 X-Y GRID RROUND RN INITIRLLY 
3-D GROSS IRN; CLOUD OF NEUTRRLLY BUOYRNT MRRKER PRRTICLES 
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FIGURE 8.4A Long Term Transport of Neutrally Buoyant Particle at Brown's Ledte with Constant Horizontal 
and Vertical Dispersion Coefficients of 50, 0, 0.0001 FT^/SEC. Two Dimensional Plot of the 
Expansion and Translation of the Horizontal Grid 
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WRPIC 2-D CONCENTRRTION CONTOURS BT VRRIOUS Z-LEVELS 


WRSTE RELEBSE SIMULRTION BT BROWNS LEDGE 
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FIGURE 8.4B Long Term Transport o£ Neutrally Buoyant Particle at :^oWn's Ledge with Constant Horizontal 
and vertical Dispersion Coefficient of 50, 0.0001 FT /SEC, . T — 0.0 Secs, Horizontal Con- 
centration Contours 
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PIGUBE 8.4C Long Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Constant Horizontal 
and vertical Dispersion Coef f icient ' of 50, 0.0001 FT /SEC. T = 3478 Secs. Horizontal Con- 
centration Contours 
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WPPIC 2-0 CONCENTRBT ION CONTOURS BT VBRIOUS Z-LEVELS 
WBSTE RELEfiSE SIMULBTION BT BROWNS LEDGE 
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8.4D Long Term Transport of Neutrally Buoyant Particle at Br,own's Ledge with Constant Horizontal 
and Vertical Dispersion Coefficient of 50, 0.0001 FT^/SEC, T = 6707 Sec. Horizontal Con- 
centration Contours 
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FIGURE 8. 4E Long Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Constant Horizontal 
and Vertical Dispersion Coefficient of 50, 0.0001 PT^/SEC. T = 12420 Sec. Horizintal Con- 
centration Contours 
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figure: 8.4F Long Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Constant Horizontal 
and Vertical Dispersion Coefficient of 50, 0.0001 FT /SEC. T =15028 Sec. Horizontal Con- 
centration Contours 
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WPP1C 2-D CONCENTRBTION CONTOURS BT VRRLOUS Z'-LEVELS 
WBSTE RELEBSE SINULRTION RT BROWNS LEDGE 
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FIGURE 8.4G Long Term Transport of Neutrally Buoyant Particle at„Brown's Ledge with Contant Horizontal 
and Vertical Dispersion Coefficient of 50, 0.0001 FT /SEC. T = 18506 Sec. Horizontal Con- 
centration Contours 




153 






COUPLING WRPIC TO 2-D VERTICfiLLY RVG. HYDRODYNRMIC TIDRL MODEL 
TRRNSLRTLON RND EXPRNSION OF 15X15 X-Y GRID RROUND RN INITIRLLY 
3-D GROSS LRN CLOUD OF NEUTRRLLY BUOYRNT MRRKER PRRTICLES 
WITH CONSTRNT X-Y ,Z DIFF. COEFFICIENTS OF 100.0 rO. 000 1 FT2/SEC 
GRID MOVEMENTS RT T^O j1 j3 iS ,fi ,7 HOURS RFTER SLRCK TIDE 
X-Y DIFF. COEFF. OF 100 FT2/SEC. WITH RN ENE NET WIND DRIFT 



FIGURE 8.5A Long Terra Transport of Neutrally Buoyant Particles at Brown's Ledte with Constant Horizontal 
and Vertical Dispersion Coefficients of lOG.O, 0.001 FT^/3EC, Two Dimensional Plot of Ej^an- 
sion and Translation of Horizontal Grid 
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WPPIC -Z-D CONCENTRRTION CONTOURS RT VRRIOUS Z-LEVEL3 
NEUTR-BLLY BUGYRNT WROTE RELERSE SIMULRTION RT BROWNS LEDCE 
TIDRL DISPERSION WITH R SUPERIMPOSED NET WIND DRIFT 
CONSTRNT H0RI20NTRL DIFFUSION COEFTICIENTS OF 100.0 FT2/SEC 
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FIGURE 8. SB Long Term Transport of Neutrally Buoyant Particles at Brown's Ledge with constant Horizontal 
and Vertical Dispersion Coefficients of 100.0, 0.001 FT /SEC. Two Dimensional Plot of Expan 
sion and Translation of Horizontal Grid. T = 0.0 Sec., Concentration Contours 
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FIGURE 8. 'SC Long Term Transport of Neutrally Buoyant Particles at Brown's Ledge with Constant Horizontal 
and vertical Dispersion Coefficients of 100.0, 0.001 FT'^/SEC. Two Dimensional Plot of Expan- 
sion, and Translation of Horizontal Grid. T — 3786 Sec., Concentration Contours 






FIGURE 


wPPIC •:-DC:ONCENTRBT I ON CONTOURS RT VRRIOUS Z-LEVEL3 
NEUTRRLLY BUOYRNT WROTE RELERSE SIMULRTION RT BROWNS LEDGE 
TIDRL D ISPERS IQN WITH B SUPER LNPOSED NET WIND DRIFT ■ 
CONSTRNTHQRIZONTRL DIFFUSION COEFFICIENTS OF 100.0 FT2/SEC 

TTNE= lOTTl. 2=^ I. TlNE= 10741. 2- 2. 

FT/ GRID UNIT 
953.7 
^Y 953.7 

WRPIC GRID ORIGIN 
LRTITUDE 41.349 
LONGITUDE 71.075 

7.U0 7.50 

X X 

CONTOUR MINIMUM 1.0 


>- 



>- 


o 
in 

rvir 



TIME= 1074 1. Z- 3. TIME== 1074 1. 2= 4. INCREMENT 




X X 


250. 0 


8.5D Long Term Transport of Neutrally Buoyant Particles at Brown's Ledge with Constant Horizontal 
and Vertical Dispersion Coefficients of 100.0, 0.001 FT /SEC. Two Dimensional Plot of Expan 
sion and Translation of Horizontal Grid, T = 10741 Sec,, Concentration Contours 
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FIGURE S.6A Long Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Scale Depend ent 2 Horx 
zontal Diffusion Dispersion Coefficients, and Vertical Dispersion Coefficient =0.01 FT /SEC 
Two Dimensional Plot of Expansion and Translation of Horizontal Grid, 
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FIGOPS 8.6B 


Long 'Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Scale Dependent Hori- 
zontal Diffusion Dispersion Coefficients, and Vertical Dispersion Coefficient = 0.01 FT'^/SEC. 
Concentration Contours at T = 0.0, 433.0, 1054.0, 2047.0 Secs 








WFPLC r-P CONCt:NTRRTION CONTOURS IN RN X-Y PLRNE 

THROUGH THE MRXIMUM CONCENTRRTION LEVEL , RT R GIVEN TIME STEP 


TIDRL TRRNSPORT OF RN INITIRLLY THREE DIMENSIONRL 
CRUSSIRN CLOUD OF NEUTRRLLY BUOYRNT NRRKER PRRTICLES 

WITH HORIZONTRL SCRLE DEPENDENT DIFFUSION COEFFICIENTSBQ 

FT/GRID UNIT 


TIME= Z- 3v 




48S. 

?G2. 


3 

a 



488. 

?62. 


3 

3 


WfiPlC GRID ORIGIN 
LRTITUDE "liim 


?l. 088 

LONGITUDE 71.084 


TIME= 6394. Z= 3. 


CONTOUR MINIMUM 1.0 
INCREMENT 100.0 


FIGURE 8.6C Long Term Transport of Neutrally Buoyant Particle at Brown's Ledge with Scale Dependent Hori- 
zontal Diffusion Dispersion Coefficient, and Vertical Dispersion Coefficient = 0.01 FT /SEC. 
Concentration Contours at T = 3662, 6394 Secs. 





COUPLING NfiPlC TO VERTlCfiLLY RVG. HYDRODYNRMIC TIDRL MODEL 
TRRNSLRTION RND EXPPtNSION OF 15X15 X-Y GRID RROUND RN INITXRLLY 
3-D 6RUSS1RN CLOUD WITH R SINGLE PRRTICLE SETTLING RRTE 
WITH CONSTRNT X-Y ,Z DIFF. COEETICIENTS OF 100. .0.0001 FT2/SEC 


GRID MOVEMENT 


INITIRL X-Y TCTRL GRID SPRN OF 3000 FT. 


HOURS RFTER TIME 1^=0 . BT SLRCK TIDE 


f Sedime 
Vertical 
onal Plo' 


xcle wxth a Single Particle Settling Ra 
n Coefficients of 100.0, 0,0001 FT^/sec 
sion and Translation of Horizontal Grid 










FIGURE 


WRPIC 2-D CONCENTRRTION CONTOURS IN RN X-Y PLRNE 

THROUGH THE RRXIMUM CDNCENTRRT ION LEVEL . RT R GIVEN TIME STEP 

WRSTE RELERSE SIMULRTION RT BROWNS LEDGE 

S-D GROSS IRN CLOUD WITH R SINGLE PRRTICLE SETTLING RRTE 
PBRTICLE SETTLING VELOCITY = 0.003 FT/SEC 

FT/ GRID UNIT 

5 200.0 


TIME= 


0. Z= 3. 


TIME*= 3 08. Z= 3. 



X 



X 


$ 


250, 0 
312. 5 
390. 8 

200 . 0 
250. 0 
312.5 
390. 8 


WRPIC GRID ORIGIN 


LRTITUDE 


LONGITUDE 


41. 321 
41. 322 
41. 324 
41. 325 


?1. 092 
71. 091 
71. 090 
71,091 


TIME*= 805. Z- 3. 


TIME- 1302. Z- 3. 



X 



X 


CONTOUR MINIMUM 1.0 
INCREMENT 100.0 


cn 

lO 


8.7B Long Term Transport of Sedimentary Particle with a Single Particle Settling Rate, with Con- 
stant Horizontal and Vertical Dispersion Coefficients of 100.0, 0.0001 FT /SEC, W = 0.003 
FT/SEC, Concentration Contours at T = 0.0, 308, 805, 1302 Secs ® 






WPPIC 2-0 CQNCENTRRTION CONTOURS IN RN X-Y PLRNE 

THROUGH THE MRXIMUM CONCENTRRTION LEVEL , RT R GIVEN TIME STEP 

WRSTE RELERSE SIMULRTION RT BROWNS LEDGE 

3-D GRUSSIRN CLOUD WITH R SINGLE PRRTICLE SETTLING RRTE 
PRRTICLE SETTLING VELOCITY = 0.003 FT/SEC 

FT/ GRID UNIT 


TIME:= 22 9B. z= 4. 


TIME- 4034. Z^4. 



X 



5 


GIO. 4 
7G2. 9 
953. ? 

488. 3 
810. 4 
?62. 9 
953. 7 


WRPIC GRID ORIGIN 


X 


LRTITUDE 


LONGITUDE 


41. 328 
41. 332 
41. 343 
41. 351 


71.090 
71. 083 
71. 079 
71.068 



FIGURE 8. 7C Long Term Transport of Sedimentary Particle with a Single Particle Settling -Rate, with Con 
stant Horizontal and Vertical Dispersion Coefficients of 100,0, 0.0001 FT^ /SEC, w_ = 0.003 
FT/SEC, Concentration Contours at T = 229S, 4304, 8007, 11362 Secs ® 








FIGURE 


wRPlC 2-D CONCENTRRTION CONTOURS IN RN X-Y PLRNE 

THROUGH THE HRXIMUN CONCENTRRT ION LEVEL i RT R GIVEN TINE STEP 

WRSTE RELERSE SIMULRTION RT BROWNS LEDGE 

3-D CRUSSIRN CLOUD WITH R SINGLE PRRTICLE SETTLING RRTE 
PRRTICLE SETTLING VELOCITY »= 0. 003 FT/SEC 


FT/ GRID UNIT 

TINE- H3R4. Z=10. TINE- 18630. Z- 3. ^ 953.? 

353.? 

1192.1 

5 953.? 

Y 953,7 
953,7 

WRPIC GRID ORIGIN 

4 1. 353 

LRTITUDE 41.333 
41. 352 


71. 059 

LONGITUDE 71.049 
71.033 



TINE- 23285. Z=10. 



X 


CONTOUR NININUN i. 0 
INCRENENT 100.0 


,7D Long Term Transport of Sedimentary Particle with a Single Particle Settling ’Rate, with Con- 
stant Horizontal and Vertical Dispersion Coefficients of 100.0, 0.0001 FT^/SEC, W = 0.003 
FT/SEC. Concentration Contours at T = 14964, 18690, 23285 Sec ® 





DEPTH-FT. DEPTH-FT 


165 . 


WflPIC SIMULATION OF SEDIMENT TRANSPORT AT BROWNS LEDGE 
FALLOUT FROM AN INSTANTANEOUS GAUSSIAN RELEASE 
PARTICLE SETTLING VELOCITY^ 0.0029 FT/SEC 
X-Z AND Y-Z PLANAR VIEWS OF PARTICLE POSITIONS 



X-GRID UNITS 



SETTLING VELOCITIES FT/SEC 
0.0030 


LOCATION OF CENTER OF NASS 
LATITUpEs= m,3l5 
LONGITUDE^ 71, 035 

X-Y GRID LENGTHS FT, 

DXs= 200. 

DY:= 200. 

TIME== 0.00 HR. 


Y GRID UNITS 


figure 8.7E Long Tenn Transport of Sedimentary Particle v?ith a Single 

Particle Settling Rate, with Constant Horizontal and Vertical 
Dispersion Coefficients of 100.0, 0.0001 PT“ /SEC, = 0,003 
FT /SEC? X-Y, Y-Z Planar View of Particle Cloud, T =0.0 Hours 




DEPTH-FT. DEPTH-FT. 

- 30 . 00 - 60 . 00 - 30 . 00 - 30 . 00 - 60 . 00 - 30 . 00 


166 


WaPIC SIMULATION OF SEDIMENT TRANSPORT AT BROWNS LEDGE 
FALLOUT FROM AN INSTANTANEOUS GAUSSIAN RELEASE 
particle settling VELOCITY^ 0, 0029 FT/SEC 
X'-Z AND Y-Z PLANAR VIEWS OF PARTICLE POSITIONS 


1 1 3 X 



J 1 ( L 


3, QO 6. 00 3. 00 12. 00 

X-GRID UNITS 

—j — ^ — , j , — 


t 



4 __i 4- 1. 


3. 00 6, 00 3. 00 12. 00 


SETTLING VELOCITIES FT/SEC 
0. 0Q3O 


LOCATION OF CENTER OF NRSS 
LPTITUDE=^ 41. 324 
LONSITUDE?= 71.030 


X-Y GRID LENGTHS FT, 
DX*= 763. 

DY- 763. 

TIME^ 2.22 HR. 


Y-GRID UNITS 
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X-Z, Y-Z Planar views of Particle Cloud, T = 1.81 Hours 




183 . 


WBPIC SIMULRT I ON OF SEDIMENT TRANSPORT RT BROWNS LEDGE 

FALLOUT FROM AN INSTANTANEOUS GAUSSIAN RELEASE 

FOR FIVE DIFFERENT PARTICLE DIAMETERS 

X-Z AND Y-Z PLANAR VIEWS OF PARTICLE POSITIONS 



X-GRID UNITS 


SETTLING VELOCITIES FT/SEC 
0. 0030 
0 . 0011 
0. 0075 
0,016? 

0, 0256 


location of center of mass 

LATITUDE^ 41,323 
LONGITUDEr= 71,031 

X-Y GRID LENGTHS FT, 

PX- 610, 



DYp= 610, 

lim- 2.16 HR,- 


Ie?SLa'^l of a sadimsaf Cloud with Dlffsror 

Settlxng Velocitxes. wxth Horxzoi'ital and Vertical Diffusion Goefflri 

0.008, 0.017, 0.026 

X-Z, Y-3 Planar Vxews of Partxcle Cloud, T = 2.16 Hours 







184 


WBPIC SIMULATION OF SEOIMENT FALLOUT AT BROWNS LEDGE 
FOR A SINGULAR PARTICLE SETTLING VELOCITY 
THE X»Y»Z POSITIONS OF CENTER OF MASS, AND % OF MASS 
REMAINING IN THE PARTICLE CLOUD AS A FUNCTION OF TIME 



TIME-HOURS 



FiGrjRE 8.9B Long Ter/n Transport of a Sodimont cloud with Five Lifferont 
Settling VelocltiEs, with Horizontal and Vertical rjif fusion Coefficients of 
100.0, 0.0001 FT /SBC, W = 0.002, 0.003, 0.008, 0.017, 0.026 FT/SBC. X-Y-Z 
Positions of the Center erf Mass and Percent Mass Remaininfj in the Particle 


3 CD 



185 






















WRPIC 2-r CONCTNTRRTION CONTOURS IN RN X-Y PLRNE 

THROUGH THE- NRXIMUM CONCENTRRTION LEVEL i RT R GIVEN TIME STEP 

X »Y »Z DIFFUSION COEFFICIENTS =100. jIOO. jQ. 01 FT27SEC 

3-D ORUSSIRN CLOUD WITH 5 PRRTICLE SETTLING RRTES 
SETTLING VELOCITIES = 0. 001 »0. 003 ,0. 000 fO. 01? »0. 026 




FT/ GRID UNIT 
^ 200 , 0 ’ 

X 312 . 5 

390. e 
4 83. 3 

^ 200.0 

Y 312.3 


WRPIC GRID ORIGIN 

LRTITUDE thill 
41. 3 30 
41. 333 

I rwir.T’riirtr* 71.032 
LONGITUDE 71.091 
71.090 
71.038 



FIGDRE 8.10B Long Term Transport of a Sediment Cloud with Five Different Settling Velocities, with 

Horizontal and Vertical Diffusion Coefficients of 100.0, 0.01 FT /SEC, W = 0*001, 0,003 
0.003, 0.017, 0.025 FT/SEC. Concentration Contours T = 0.0, 1054, 3910 fee 
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WIC SIMULATION OF ""SL'CllMENT TRANSPORT AT CRONNS LPDCE 

FALLOUT FROM AN INSTANTANEOUS GAUSSIAN RELEASE 

FOR FIVE different PARTICLE DIAMETERS 
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FIGURE 8.10D Long Terra Transport of a Sediment Cloud with Five Different 
Settling Velocities, with Horizontal and Vertical Diffrugion Coeff Lcients of 
100.0, 0.01 FT^/SFG, W >-0.001, 0.003, 0.008, 0.017, 0.026 FT /SEC. X-Y, 
Y-Z Planar Views of Particle Cloud, T -0.0 Hour 
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FIGURE 8.10E Long Term Transport of a Sediment Cloud with Five Different 
Settling Velocities, with Horizontal and Vertical Diffusion Coefficients of 
100.0, 0.01 FT^/SEC, W - 0.001, 0.003, 0.008, 0.017, 0.026 FT/SEC. X-Z, 
Y-Z Planar Views of Particle Cloud, T = 1,09 Hour 
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FIGURE 8. lOF Long Term Transport of a Sediment GlonS with Five Different 
Settling Velocities, with Horizontal and Vertioal Diffusion Coefficients 
of 100.0, 0.01 FTVSEC, W = 0.001, 0.003, 0.008, 0.017, 0,026 FT/SEC. X-Z 
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WRPIC SinULRT I ON OF SEDIMENT FRLLOUT RT BROWNS LEDGE 
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FIGUlffi 8,10G Long Terra Transport; of a Sediment Cloud with Five Different 
settling Velocities, with Horizontal and Vertical Diffusion Coefficients 
of 100.0, 0.01 PT'^/SEC, W = 0.001, 0.003, 0.008, 0.017, 0.026 FT/SEC. X-Y-Z 
Positions of the Center of Mass and Percent Mass Remciining in the Particle 
Cloud as a FunGtipn of ; - 
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IX. SOMHARY, CONCLUSIONS, RECOMNEND ATIONS 
WAPIC, a three dimensional Lagrangian marker particle 
in Eulerian finite difference cell technique for simulating 
the incompressible mass transport equation, has been 
developed. The method was verified against five analytic 
solutions and subsequently applied to the prediction of 
long term tidal transport of suspended sediment from an 
instantaneous barge disposal of dredged sediment waste at 
Bro wn ’ s Ledge. 

A. Verification Studies 

The procedure was found, when utilizing a horizontal 

standard deviation to grid spacing ratio of 2.0, to 

accurately predict to within an absolute value of 6 % , the 
analytic concentration values for three dimensional 

unbounded turbulent transport of an initially Gaussian 
release in uniform steady unidirectional flow and steady 
uniformly sheared unidirectional flow. 

WAPIC was applied in the simulation of long term 
transport of continuous Gaussian releases in a steady 
unidirectional turbulent flow field. comparison with 
analytic concentration values were poor because of 

inadequate numbers of marker particles to represent the 
initial Gaussian particle distributions. In general, 
WAPIC, in its present form, is impractical for simulations 
of continuous waste releases, because of the excessive 
amounts of computer or auxiliary storage required for the 
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large number ;of marker particles necessary for adequate 
spatial represeritation of the succession of Qaussian 
releases* An accurate representation of the diffusive 
behavior of the continuous waste release may require on the 
order of 30 to 50 thousand marker particles in order to 
perform a two hour real time simlation, 

WAPIC was verified against an analytic solution for 
turbulent transport of an instantaneous release in uniform 
plug channel flow. Predicted concentration values were 
within 5 X of the analytic values in the center of the 
channel. This simulation demonstrated the utility of 
WAPIC*S reflection boundary conditions and procedures for 
selective translation and expansion of the Eulerian finite 
difference grid. 

WAPIC successfully simulated Dobbins' (73) one 
dimensional mass transport solution "for fallout of an 
initially vertically uniform sediment distribution with a 
single vertical settling velocity in a steady turbulent 
channel flow field. The significance of this simulation 
was the verification of WAPIC'S ability to accurately 
predict the vertical concentration and fallout rates of 
neqatively buoyant particles in a turbulent flow field. 
Vertical concentration values and particle deposition rates 
were all within 5 % accuracy for five different particle 
settling velocity cases. The Courant dynamic stability or 
time step restriction which is utilised in Eulerian finite 
difference schemes, was found to be a valid accuracy 
requirement for the WAPIC Eulerian finite difference - 
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Laqrangian marker particle algorithm. 
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B. Simulations At Brown’s Ledge 
Seven simulations of long term transport of suspended 
sediment resulting from an instantaneous disposal of 
dredged waste at Brown’s Ledge in Rhode Island Sound were 
performed, a two dimensional vertically averaged tidal 
hydrodynamic model was successfully coupled through the use 
of a velocity processor to WAPIC. For all seven neutrally 
and non-neutrally buoyant suspended sediment simulations, 
the initial plume of hypothetical maximum concentrations of 
1500 to 3000 mg/liter dropped to ambient background 
concentration levels of 1-5 mg/liter in less than two tidal 
cycles in approximately 60 to 120 feet of water. 

The center of mass of the turbid cloud for all four 
negatively buoyant simulations was "conservatively" placed 
relatively close to the surface at a depth of twenty feet* 

In reality, the center of mass of the turbulent cloud would 
probably be located much closer to the bottom. The 
suspended sediment cloud would be produced by turbulent 
entrainment or mixing with the water column as the initial 
release descends and by the vertical upward momentum surge 
of the spoil cloud following the impact at the bottom (77) . 

For the given release area shown, provided the vertical 
turbulent diffusion coefficients and free fall velocities 
are representative, the entrained cloud of sediment Would 
settle to the bottom or diffuse to background levels in a 
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much shorter time span than shown in these simulations. 
The particle clouds were also conveniently placed near the 


surface for visual demonstration of the vertical descent of 
the clouds to the hottcm. 

The most important aspect of these simulations is that 
the long term horizontal and vertical transport of a 
spectrum of particle settling velocities can now be 
modeled, provided the ambient current field, horizontal and 
vertical turbulence levels, and physical characteristics of 
the suspended sediment are accurately known. 


C. Fecommendations 

HAPIC, in its present form, may simulate the long term 
transport of an instantaneous release of sewage or 
pollutants from an outfall in a tidally dominated estuarine 
environment. Of interest would be the coupling of MAPIC to 
a three dimensional finite difference advectlve velocity 
field. Gordon (36) has developed a three dimensional 
finite difference hydrodynamic model for the Providence 
River, Utilizing the strictly Eulerian (non translating 
and expanding) version of WAPIC, the two models may be 
combined for the prediction of pollutant transport. At 
present the WAPIC velocity interpolation scheme is only two 
dimensional. Should WAPIC be coupled to a three 
dimensional hydrodynamic field, and the Eulerian grid 
system expands and translates in either two or three 
dimensions, then a three dimensional velocity interpolation 


schene will be required. 

A field verification of the WAPIC mass transport 
procedure in a tidally driven environment is necessary in 
order to deteriDine the actual predictive capabilities. A 
search for a suitable field study of long term pollutant or 
sediment transport with an adequate data base of acivective 
field velocity vectors, tide heights, and suspended 
constituent levels for testing the procedure should be 
undertaken, 

WAPIC is suitable for predicting the long term 
transport of dilute suspended sediment, but does not 
address the short term aspects of the descent, impact, and 
spreading surge of the negatively buoyant sediment or 
pollutant cloud. A two phase short term hydrodynamic model 
would be of value to produce the initial sediment 
distribution for WAPIC*S long term transport simulation. A 
suitable simulation technique might be a hybrid ALE, ( 
Alternating Eulerian Lagrangian) , PIC, ( Particle In Cell) , 
incompressible hydrodynamic technique. The PIC method is 
well suited for simulating the hydrodynamics of 
incompressible two phase flow situations. 

WAPIC may be utilized for simulating the erosion and 
entrainment of sediment from the bottom boundary layer. A 
” Shields'* critical shear velocity criteria (78) for 
entrainment of sediment might be considered , The 
vertically transformed, horizontally expanding and 
translating version of WAPIC coupled to a two or three 






JHl 


laa g HB r g« j ! a» » gja ig BKa B8e aiiHi^MW^^ 


197 . 

dimensional hydrodynaniic velocity field might be used to 
study the sediment erosion and deposition patterns within 


an estuarine environment 
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Derivation of transformation identities (36) 
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Siibstituting for c, we define the derivative of A with respect to x 
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_ z - §(x>y>t) 
H (x, y, t) 


where 


H(x, y, t) = §(x, y, t) + h(x, y) 


£!1 

fiz 


[ 

6z *- 


z - § _ 1 

'~1T~ J ^ H 


bz 

bt 




5 

St 


+ 1)§ + T] h] = (T] + 1) || 


0 

0 

0 

. Y, Tl, t 

(. 

V 

V 
C 

( 




APPENDIX A 


CONTINUED 


SO 


6t 


bt 


* I i T 


Prom A. 2 for x-transformation 


6 A ^ b_A ^ 6T] bA 

ftX 6z fix 


^ = r- + i* [ ( T) + 1) ^ + Tl ^ ^ 

fix bx H ^ ' V fix . fix bTl 


Similiarly for Y-transformation 


+ i r /"n + 1 ^ + 'll ^ 

TV - ~ + H L<"f^ + 6y + ^ ^ ^ b1T 


fiy 5y H 
For z-term Transformation 


iU = -L. r 2-1-1 1 M = 1 

' bll "bz ^ H '* 


b^T " H b^I 


(A.9) 


(A.IO) 


(A.ll) 


{A. 12) 


(A. 13) 
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APPENDIX B 

VERTICAL TRANSFORMATION OP MASS TRANSPORT EQUATION 

The transformation of the cartesian three dimensional mass transport 
equation may be accomplished using the transformation identities A. 9, 

11, 12, 14 after Gordon (36). For convenience, the 3-D mass transport 
equation will be divided into seven (7) terms, 


L° + 

bt 

n2.° + 

bx 

V ^ + 

by 

bz 

- -L_ (K 

bx X 

(a) 

(b) 

(C) 

(d) 

(e) 


bc 


) - 


-^(K 

by y bz 

(f) 


^(K |^)=0 

bz Z 


(g) 


each term being transformed separately. 


(B.l) 


Li = ^ IRT] + 1) ^ 


6 tJ bT] 

= U[ — _ i [(1 + 7]] ^ + T] ^ -j Li 3 
bx ^ fix H fix fi^ J b^ 


(B.2) 


(B.3) 


V Li = v[ ^ 
by fiy H 


- i [ (1 + Ti) ] LL + Ti ~ ^ ] 

H ^ ^ fiy fiy ■* bTl 


(B.4) 


Consider the total time derivative 


d 

dt 


= h... 

bt 


+ u + V ^ 


b^ 


by 


+ w 


bz 


(B.5) 


Transforming this 


d 

'^t 


= + u 

fit 


+ V i -I- - i -( Li 

fcy H Vff 5t 


° ^ r" I > w 


Defining a dimensionless vertical velocity 


W 1 ^ 

H “ H ^ bt 


+ U Li + V Li ) 
. by - 


(B.7) 


0) 
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APPENDIX B - EQUATIONS - CONTINUED 


Thus 


W = 


CUH+(^ + u^ + V~) 

Ot 6x ' 


(B.8) 


” ^ “ [«^H + ( ^ + u ~ + V ~) 1 i. ^ 

52 Ot 5x Oy' •» H ^T\ 


(B.9) 


-2 -[k i£] 

IX X bx-* 


^;k - i ri + Til ^ ^ 1 ^ 1 

^ X>-6X H L .X ' *x J xTl 


6* ' feT] 


+ HyO. T 


^ [K ii-- 1 
Oy L y 5y ' 


-.§,.. r ^ i. rl 4. TlT 4. 41 -i b *^ 1 

6 Y ^ fiY H C fiy 6Y ^ bTli 


(B.IO) 


(B.ll) 



+ T. 


K 


[K ^ ] = -IL. ^ 

b^ z 52 bTl ^2 j^Tl 


(B.12) 


Now the transformed terms in each dimension will be rearranged in each dimen- 
sion. 


[x - terms^ - B. 3 without u term and B. 10 

bT| 


U 


fix 


-1- [k [ - in +111 T ^ 

6x '-V 6X H LI + ’ll It T 


(B.13) 


assuming incompressibility expression B.13 may be rearranged as 
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APPENDIX B - EQUATIONS - CONTINUED 


Jl. f- 

6X c 6x 


+ Ti||]a£]c 


(B.14) 


5X T 


where 


(B.15) 


u = u - ~ ^ - I [1 + TH + T1 f - 1 

T G gx H. fix bM 


(B.16) 


j-y-Terms^ Similarly 


V„ = V 


r- I- M L I e^r fiY ■* bTl 


(B.17) 


~ (V C) 

6V T 


(B.18) 


[jl - Terms] 


Combining ^ terms of B. 2, 3, 4 with B.9 


Assuming 


^ , ■ b£ , hi , 2)^ , il , hi 
bt bt 5Y bx by 


< < 




with some rearrangement 


±^r\,lll hi 
5T1 L ^ 


1 ru ^ V ^ 1 ] c 
H L fiX 6Y 


(B.19) 


which is equivalent to 


fi T] ^ 


(B.20) 


ConbiTixng B.15, 17, 20, 21 


+ -T- [W^c"' 

6y ^ T 


[tu^c" = 0 


where U^, V^, and cu^ are defined by B.16, 18, and 20, 

B.22 represents the vertically transfor'ned total pseudo velocity form of 


the three dimensional mass transport equation. 
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